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I. EXECUTIVE SUMMARY 


The imminent destructive threats of Lightning on helicopters and other airborne sys- 
tems has always been a topic of great interest to this research grant. Previously, the 
lightning induced currents on the surface of the fuselage and its interior were pre- 
dicted using the finite-difference time-domain (FDTD) method as well as the NEC 
code. The limitations of both methods, as applied to lightning, were identified and 
extensively discussed in the last meeting. After a thorough investigation of the ca- 
pabilities of the FDTD, it was decided to incorporate into the numerical method a 
subcell model to accurately represent current diffusion through conducting materials 
of high conductivity and finite thickness. Because of the complexity of the model, its 
validity will be first tested for a one-dimensional FDTD problem. Although results 
are not available yet, the theory and formulation of the subcell model are presented 
and discussed here to a certain degree. 

Besides lightning induced currents in the interior of an aircraft, penetration of 
electromagnetic fields through apertures (e.g., windows and cracks) could also be 
devastating for the navigation equipment, electronics, and communications systems 
in general. The main focus of this study is understanding and quantifying field 
penetration through apertures. The simulation is done using the FDTD method 
and the predictions are compared with measurements and moment method solutions 
obtained from the NASA Langley Research Center. 

Cavity-backed slot (CBS) antennas or slot antennas in general have many appli- 
cations in aircraft-satellite type of communications. These can be flushmounted on 
the surface of the fuselage and, therefore, they retain the aerodynamic shape of the 
aircraft. In the past, input impedance and radiation patterns of CBS antennas were 
computed using a hybrid FEM/MoM code. The analysis is now extended to coupling 
between two identical slot antennas mounted on the same structure. The predictions 
are performed using both the hybrid FEM/MoM and the FDTD NEWS code. The 
results are compared with each other as well as with measurements performed in the 
ElectroMagnetic Anechoic Chamber (EMAC) of ASU. In addition to self and mu- 
tual impedances versus frequency, the comparisons include mutual coupling Si 2 as a 
function of distance for various slot orientations. 

The FDTD NEWS code was upgraded from second-order accurate in time and 
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space, e.g. FDTD(2,2), to second-order accurate in time and fourth-order accurate 
in space, e.g. FDTD(2,4). It was shown in the previous report that the higher-order 
FDTD scheme is more accurate for the same discretization size and exhibits smaller 
dispersion errors. This was illustrated for a one-dimensional wave propagation prob- 
lem. In this report, an in-house developed three-dimensional version of FDTD(2,4) 
is used to predict radiation patterns of electrically large helicopter problems. Specif- 
ically, the radiation patterns of a monopole on the NASA scale helicopter model at 
9.18 GHz are predicted and compared with the FDTD(2,2) NEWS code and mea- 
surements. 

A new iterative algorithm is formulated in conjunction with the hybrid FEM/MoM 
approach to solve effectively for the coupling parameters of multiple cavity-backed slot 
antennas on the surface of a platform. Each radiating element is solved independently 
to find the governing field distribution in the aperture of the cavity. The interaction 
among these cavities is accounted for through an iterative procedure which continu- 
ously updates the fields in each finite element domain until convergence is achieved. 
The validity of the method is tested by calculating the mutual coupling between two 
similar cavity-backed slot antennas mounted on an infinite ground plane. This ap- 
proach is not only accurate but also extremely powerful and efficient both in terms 
of memory storage and solution time. 

The hybridization of FEM with high-frequency methods such as the physical optics 
(PO) has been a focus of this project for the last few years or so. In this report, the 
FEM is formulated for two-dimensional scattering and radiation problems in conjunc- 
tion with either absorbing boundary conditions or the boundary integral approach. 
The validity of the method has been tested against MoM solutions. For radiation 
problems, the finite element domain is coupled (using a one-way interaction) to a 
physical optics formulation to compute scattered fields by a nearby large object. 

Spectral methods are applied for the first time to the solution of one-dimensional 
problems. The great potential of these methods in the area of electromagnetics is 
demonstrated through numerical experimentation. Conclusions are obtained by com- 
paring results with the standard FDTD algorithm. Furthermore, the PML concept 
has been applied to spectral methods to simulate open-space propagation problems. 
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II. INTRODUCTION 


The major research topics addressed in this report are the ones which have been 
identified and recommended by the Advisory Task Force of the Advanced Helicopter 
Electromagnetics (AHE) program. Although some of these topics provide a contin- 
uation of previous work, emphasis was placed on research areas that are of special 
interest to the program. Some of these subjects were brought to our attention in the 
Annual Conference held at the Boeing plant in Philadelphia, PA, on May 26-27, 1999. 
The main topics of focus in this report are the following: 

• Lightning and its numerical modeling 

• Field penetration through apertures 

• Analysis of cavity-backed slot antennas 

• Mutual coupling of cavity-backed slot antennas 

• Higher-order FDTD schemes 

• Hybrid methods 

• Spectral methods 

Each of the above topics is addressed in sequence in this document, reporting on 
progress already made primarily during the period after the May 1999 meeting and 
outlining future work to be accomplished. 

The topic of lightning is now being studied in a more systematic way. Two different 
numerical methods have been used to predict the diffusion of currents in the interior 
of a conducting enclosure of finite skin thickness and conductivity. One of these two 
methods was the NEC code which models the structure using a wire-grid model. The 
thickness of the walls was taken into account implicitly through the radius of the 
wires. Although the results predicted by NEC agreed well with results obtained from 
an independent study performed by Dr. Jack Nachamkin at Boeing in Philadelphia, 
the NEC analysis does not account for the skin depth of the wall, which is known to 
be a function of frequency and conductivity. A more accurate approach to predict the 
diffusion of currents in the interior of a conducting structure is to use a more robust 
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technique such as the FDTD method. The FDTD method is accurate provided the 
discretization inside the thickness of the walls is adequate enough. Unfortunately, it 
was determined in previous reports that the cell size of the grid inside the conductor 
has to be extremely small to achieve accuracy in the solution. In addition, the larger 
the conductivity of the wall, the smaller the cell size should be. This will increase the 
computational domain substantially, therefore limiting the use of the FDTD to only 
small problems. 

In this report, the FDTD scheme is augmented by introducing a subcell model 
to accurately represent field and current variation inside the thickness of the walls. 
This subcell model assumes a one-dimensional variation of the fields of the slab in a 
direction normal to the wall. This assumption is valid provided, first, the thickness 
of the wall is much smaller than its radius of curvature and its transversal direction 
and, second, the conductivity of the material is sufficiently high so that the waves 
propagate in a direction normal to the wall. Because of its complexity, this model 
will be first implemented for a one-dimensional problem. At this stage, numerical 
results are not available; thus, emphasis will be directed toward the formulation and 
the theory behind the subcell model. 

Penetration of field intensity into the the interior of an aircraft through apertures 
is an EMI issue with great relevance to the AHE project. When an aircraft is fly- 
ing near high-power transmitting antennas, the external threat might couple through 
apertures with sensitive control equipment thereby resulting in severe malfunction- 
ing. The ability to predict the amount of power that couples through an aperture 
at a given frequency is of great interest to us. In this report, the FDTD method is 
used to compute field penetration inside a moderately sized cavity with three aper- 
ture configurations: a single large aperture, two large apertures, and a single small 
aperture. The object is excited with a normally incident Gaussian pulse with spec- 
trum up to 1 GHz. The electric field shielding, which is the relation between the 
field in the interior of the box relative to the field outside the box, is plotted for all 
three cavity configurations as a function of frequency. The FDTD predictions are 
compared with numerical data obtained using the method of moments. Although at 
first this problem might seem trivial for a method such as the FDTD, there are for 
sure numerical difficulties associated with the highly resonant nature of the cavity as 
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well as the electrically small size of the aperture. These problems can be alleviated by 
introducing a loss mechanism into the computational domain to accelerate the decay 
of the late-time fields. 

Antenna technology as applied to helicopters was always of special interest to 
the AHE project. During the last few years, emphasis was placed on the numerical 
analysis of cavity-backed slot (CBS) antennas which can be conveniently flushmounted 
on the surface of a fuselage thereby retaining the aerodynamic profile of the aircraft. 
Prediction of input impedance and radiation patterns of a single CBS antenna on 
an infinite ground plane was achieved through a hybridization of the FEM and the 
spectral/spatial domain method of moments. This hybridization offers numerous 
advantages since the free-space region does not need to be descritized; this saves 
memory space and speeds up computational time. The FDTD is also used in this 
report to predict the input impedance of a single CBS antenna. Although this problem 
may seem trivial, it is actually difficult since the time-domain fields propagating 
inside the cavity decay with a very slow rate. As is discussed in this report, the 
source modeling has been augmented with a source resistance which helps dissipate 
the time-domain fields inside the highly resonant structure. 

Both the hybrid FEM/MoM method and the FDTD method have then been ap- 
plied to compute the mutual coupling of two identical CBS antennas mounted onto 
a ground plane. For such a problem, the FDTD has a severe drawback compared to 
the hybrid approach. Using the FDTD, the free-space exterior-to-the-cavities region 
needs to be discretized and terminated properly with some type of absorbing bound- 
ary conditions. This approach creates large computational domains which, in turn, 
generate large errors primarily due to dispersion. The hybrid FEM/MoM approach, 
on the other hand, does not discretize the exterior region; thus, placing the two CBS 
antennas farther away presents no significant challenges to the method. Predictions 
on mutual coupling using the hybrid FEM/MoM approach are shown versus frequency 
and distance and compared with measurements performed in the EMAC facility of 
ASU. 

This hybrid FEM/MoM approach was then improved to further speed up the 
computational time as well as to further reduce memory requirements. Specifically, 
instead of treating the two CBS antennas as a tightly coupled system which can be 
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solved directly using the hybrid method, they are rather perceived as two non-coupled 
antennas placed in close proximity to each other. For mutual coupling computations, 
one antenna is always excited using a current or voltage source, whereas the second 
antenna is left as an open circuit. The field distribution for the transmitting antenna 
is obtained in the absence of the second antenna. Through a coupling matrix, which 
basically represents the governing transfer function, the induced field in the aperture 
of the second antenna is computed. This field acts as an excitation to the second 
antenna problem which can be solved to obtain the field distribution everywhere 
including the cavity. Using this field, a first-order approximation to the mutual 
impedance may be computed. The same coupling matrix can be used to compute the 
induced by the second antenna field in the aperture of the first antenna. The iterative 
process repeats until convergence is achieved. Results based on this algorithm are 
presented in this report for the coupling of two identical CBS antennas placed at 
some distance apart. The predictions are compared with data obtained using the 
hybrid FEM/MoM direct approach. Results of mutual coupling as a function of 
distance are also included in the report. 

Higher-order techniques in the context of the FDTD method have been one of our 
research topics during the last year. In the last couple of reports, we have introduced 
the concept of higher-order FDTD schemes. Specifically, a second-order accurate 
scheme in time and fourth-order accurate in space [FDTD(2,4)] was introduced for 
the solution of one-dimensional wave propagation problems. The results obtained 
using FDTD(‘2,4) were compared with results from a standard FDTD(2,2) scheme. 
It was shown in a previous report that the higher-order scheme exhibits significantly 
lower dispersion error for the same discretization. As a result of this observation, 
the higher-order scheme can be used to more accurately solve for electrically large 
problems, which are of special interest to the AHE project. Although there may seem 
at first that the implementation of FDTD(2,4) into a generic three-dimensional code 
is straightforward, there are few challenging problems that still need to be overcome. 
For example, for a fourth-order accurate scheme in space, the accurate enforcement 
of Dirichlet boundary conditions on PEC walls and around discontinuities is not easy. 

The FDTD(2,4) has recently been implemented into a three-dimensional code to 
more accurately solve radiation problems at higher frequencies. For the truncation 
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of the computational domain, the anisotropic PML was used. This new higher-order 
code is very similar to the NEWS FDTD(2,2) code with the exception that scattering 
has not yet been implemented. This is the first version of FDTD(2,4) which is still 
under validation. As part of this validation process, the radiation patterns of a 
monopole mounted on the NASA helicopter are computed at 9.18 GHz. Numerical 
results are compared with the FDTD(2,2) (NEWS code) method and measurements. 

For electrically large radiation problems, it is almost impossible though to obtain 
an accurate solution within a reasonable amount of time. One solution to electrically 
large problems is the implementation of a robust hybrid method. In this report, we 
have already introduced a hybrid FEM/MoM approach which computes mutual cou- 
pling through an iterative mechanism. The same ideas in now under investigation for 
the hybridization of FEM with PO to solve antenna radiation problems for helicopter 
communications. In other words, the fields of the transmitting antenna will be solved 
at the first iteration using the FEM. The PO surface currents on the surface of the 
fuselage will be computed based on the radiated fields by the antenna. These surface 
currents produce scattered fields that are affecting the first-iteration field intensity in 
the vicinity of the antenna. This perturbation field will be considered at the second 
iteration to correct the governing fields in the antenna domain. The iteration contin- 
ues on until convergence is achieved. To implement this idea, it was decided to begin 
with a two-dimensional finite element code which uses either ABC’s or BI methods 
for the truncation of the computational domain. This code has been written and 
validated against other numerical solutions in this report. A first-order hybridization 
with the PO is also shown here. This means that, at this stage, there is only one-way 
interaction between the antenna and the object; i.e., the fields scattered by the object 
are not allowed yet to change the field distribution in the vicinity of the antenna. This 
feature will be implemented in the next AHE report. 

An alternative to hybrid methods approach to solving electrically large computa- 
tional problems is the use of spectral methods. Spectral methods are very powerful 
due to their impressive accuracy, exponential convergence, and negligible dispersion 
and dissipation. These are all desirable features for the solution of radiation prob- 
lems involving helicopters. In the previous report, we presented and discussed the 
basic formulation of spectral methods. Some of the feature of spectral methods were 
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demonstrated through numerical experiments and comparison with finite difference 
methods. In this report, the spectral methods are applied for the first time to one- 
dimensional propagation problems. The accuracy and future potential of these meth- 
ods in electromagnetics is demonstrated through numerical experiments. 
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Chapter 1 

Thin Conducting Layer Subcell 
Modeling 

I. Introduction 

Problems concerning numerical accuracy may arise when someone tries to simulate, 
using FDTD, structures that contain elements such as sheets, slots and wires which are 
too fine to be resolved by the smallest affordable cell size. The problem becomes more 
serious when the material is of high conductivity which must be distinguished from a 
perfect conductor. In this section, we present a subcell modeling for thin conducting 
sheets [1], This model assumes a one dimensional variation of the fields within the 
slab in a direction normal to the sheet. This assumption is a good approximation 
under the following conditions: 

• The thickness of the sheet is much smaller than its radius of curvature and its 
transversal dimensions 

• The conductivity of the material is sufficiently high so that the waves propagate 
in a direction normal to the sheet. 

II. Analysis 

We start the analysis by examining the role of the second condition mentioned in 
section I.. For this purpose, the problem in Fig. 1.1 is examined. In this figure, a 
uniform plane wave is incident from free space on a conductive material half space. 
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Figure 1 . 1 : Wave incident from free space on a conductive half space 
The incident wave electric field E t is described by: 


Ei = a z E Q e- ]0o{xcos9 ‘ +ysin8 ' ] = a z E 0 e~^ f 

(i.i) 

where 

A> = Vy/HotO 

The transmitted wave is given by: 

( 1 . 2 ) 

E t = a z E t oe -T»*-7W = a z E t Q e~^ f 

(1.3) 

The coefficients 7 *, 7 y satisfy the constraint equation: 


ll + = 7 2 = (l - 

(1.4) 

Since the continuity of the tangential fields at the interface x — 0 must hold for any 
value of y, the coefficients of y at the exponentials of the incident ( 1 . 1 ), reflected and 
transmitted ( 1 . 3 ) fields must be the same; hence: 

7y — - J — j {3 o sin B i 

(1.5) 
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From (1.4) and (1.5) we obtain 7* : 


lx = \Jl 2 - (l - j^j + fio sin 2 0,- (1-6) 

For a highly conducting material (cr we), (1.6) and (1.2) give: 

lx = v5 J ' U<7 + $> sin2 


1 +j 

V2 


\jL0ju7 


\ 




(1.7) 


If <7 >• weo and /.i is of the same order as /zo then by writing 7* = a x + jfi x , we obtain: 

a x = fi x = a/— 77— (1*S) 


From Fig. 1.1, it is seen that t&nipt — fiylfix- Hence, under the same assumptions as 
above and using (1.8), (1.5) and (1.2), we obtain: 


tan 


fiy „ y/2 — sin - 
fix - \J^i iry 



(1.9) 


From (1.5) and (1.9), it is seen that not only the constant amplitude planes but also 
the constant phase planes in Fig. 1.1 are approximately parallel to the interface. This 
means that the fields within the conducting material vary only in the x-direction, 
which is normal to the interface. The same is true for the case of the conducting 
layer shown in Fig. 1.2. Hence, the electric field within the layer satisfies the one 
dimensional wave equation: 


\/ 2 E~i 2 E - 0 & 


E = A~<? x + /l+e"^ 


where A , .4 + are arbitrary constants and 7 is given by: 


7 = 




( 1 . 10 ) 


(1.11) 
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Figure 1.2: Thin conductive layer 


The square root in (1.11) is the one with the positive imaginary part so that the first 
term in (1.10) represents a wave traveling in the negative x-direction and the second 
a wave traveling in the positive. Using the Maxwell- Faraday equation combined with 
(1.10), we obtain the magnetic field: 

H = - — v xE 

= —a x x f A~e yx - A+e-'**) = H t (1.12) 

Uf.1 v ' 

This field has only components along y and z axes, and is denoted as H t . The purpose 
of the present derivation is to express the electric fields at the two faces of the slab 
shown in Fig. 1.2 in terms of the magnetic fields, thus producing a “two-port network” 
representation of the slab in Fourier transform domain. In order to achieve this, the 
electric and magnetic fields in (1.10) and (1.12) are evaluated at the two faces of the 
slab (x=0 and x=d): 

E (0) = A- + A+ (1.13) 


12 


E(d) = + 

(1.14) 

H (0) = H n = —d x x( A~-A + ) 

C Oj.1 v 7 

(1.15) 

H (d) = H t2 = —a x x (A~e yd - A + e~' id ) 

' Wfl V 7 

(1.16) 


where the subscript t indicates that the magnetic fields have only components tan- 
gential to the faces of the slab. Taking the cross product of a x with both sides of 
(1.15) gives: 


a x x H n = ~ a x x («® x (A — '4 + )) 




31 l -A 


U)fl 


i A t ~ Aj) 


(1.17) 


where Af, A t are the tangential components of A + , A respectively. The same 
procedure but for (1.16) gives: 


a x x H t2 = g- (A+e-v* - A;e^) 

Solving the system of (1.17) and (1.18), we find: 

ujfi ~^ d (o* x H n ) A a x x H a 


Af = ~r~ 

31 


e~~f d - e^ d 


(1.18) 


(1.19) 


LOfi ~ e ld ( a ® x Hn) + dx x H t 2 

1 ~ j7 e _7d - e~* d 


( 1 . 20 ) 


Substituting (1-19) and (1.20) into (1.13) and (1.14), the desired two-port network 
relationship is obtained: 


' E tl ' 


' z„ 

Z 12 


n x Hn 

E t 2 


Z 21 

Z 22 


i 

c* 

X 

1C 
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where, 


Z\\ — Z 21 — 


22 


JW/i 


1 


k tan (kd) 


( 1 . 22 ) 


Z\2 — Zvi — — 




1 


, . , I/A (1-23) 

k sm (kd) 

and k is such that 7 = jk so that the expressions in ( 1 . 22 ) and (1.23) are of the same 
form as in [ 1 ]. The unit vector n is normal to the slab and directed outwards in both 
faces: 


_ _ J — a x for x = 0 
] a x for x = d 

For good conductors ( a ^ use), k in (1.4) becomes: 


-31 = 




e ( 1 - j—j - = (1 - 3 ) 

' LOC J 




(1.24) 


(1.25) 


k 2 = —jusficr (1.26) 

The implementation of (1.21) in time domain starts by expanding the trigonometric 

functions of (1.22) and (1.23) in power series and keeping a finite number of terms: 

(in • y^°° (foj) 2n (-i) n 

_ JCJ/i COS {kcl) _ JUfl Z^n—O (2n)! 

11 k sin (kd) k y^°o (fcrf) 2 n+ 1 (-i) w _ 

v ' Z-^n - 0 (2n+l)! 


J± W 
k 2 d 


y^co (kd) 2n ( — l) n 
z-^n = 0 (2 n)! 

y^x> (kd) 2n {-l) n 
l^n-0 (2n+l)! 


y^oo (kd) 2n (-l) n 
1 (2n)! 

ad y^°° (ferf) 2 n (~i) n 

2^n=0 (2n+l)! 


_l_ 

ai KU 


(-ju/ 2 ff(i 2 )"(-l)" 

( 2 »)! 

(2n+l)! 


y-cc (i^d 2 )” 
I 2^n—0 (2n)! 

<j(/ y^oo OWtrrf 2 ) n 
Z^n=0 (2n+l)! 
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y'oo Q VT 
1 Z^n=0 (2rc)! 

ad n (i^T- 

2^n= 0 (2n+l)! 


yM QV)" 
i L-, n =0 (2 n)! 


where, 

Similarly, Z 12 is written as: 


ad y N iEEll 

Z-n=0 (2n+l)! 

u/ = ufiad 2 


Z 12 ~ — : 


712 “ crd bV) n 

2^ n =0 (2n+l)! 

Now (1.27) and (1.29) are expanded in a partial fraction sum: 


y- Cri^k jV + -Pil.fc 

jw' + ^ (ju/) 2 + E k jus' + Ffc 


1 [A v4ll,n 

Zl1 = ^2 -— 7 


v _ 1 I v' ^12, n , C\2,k ju' + D\2,k 

Zu - — j 2 . •J7F + 2 - 


( 1 . 27 ) 


(1.28) 


(1.29) 


(1.30) 


(1.31) 


<7tf \ “ jw' + B n fr[ (ju'f + Ek ju> + Fk J 

where the first sums correspond to the real roots in the denominator of (1.27) and 

(1.29) , and the second sums to the pairs of complex conjugate roots. The above 
expansion should be modified by including a polynomial if M > N. Substituting 

(1.30) and (1.31) into (1.21), we obtain: 


En = 


| ( rn l m l 

= —j j E 2^11, n + E ^11.* + E ^12,n + E X 

aa fn=l fc=l n= 1 fc=l 


12 ,/s 


( 1 . 32 ) 


■j ( m l m l 

Et2 = 1 E ^21, n + E -^21 , k + E 2^2 2, n + E 2 

a d ln=l Ar=l n=l A-=l 

where the auxiliary variables Xij, n , Ajyfc are defined as: 

2 I 11 n _ - 

An n = . . — x H t \ 

+ B n 


22, fc 


ATi.fc = 


_ Cn jk jui' -f D\\ t k 


{jw') + E k ju>' + F k 


n . x Hi 


ti 


(1.33) 


( 1 . 34 ) 


(1.35) 
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(1.36) 


X 


12, n 


^12, n 

ju' + 


n x H t2 


X 12lk 


C\2,k j^' + D\2,k - fy 

4 1 — n x H t 2 

(juf + E t ju' + Ft 


(1,37) 


and so on for the rest of them. These equations are transformed in time domain giving 
auxiliary differential equations. From (1.34) and (1.35), for example, we obtain using 
also (1.28): 


f-ivd 2 —Xu in (t) B n Xn in (t) — A\\ >n n x 




j2 ? 

-^XnAt) + ^d 2 E k -X u , k (t) + F k X n At) 


Her d 2 C k j t (h x H n {tj) + D k nx H n {t ) (1.39) 

These clifFerential equations are then discretized in time using a finite difference 
scheme and update equations are derived and combined with the usual Yees’ FDTD 
algorithm as follows: First, the magnetic fields Hn, H t 1 are advanced using the last 
values of E t 1 , E t 2 . Since the tangential magnetic fields at the faces of the slab are 
not included in the normal FDTD scheme, they can be approximated by the nearby 
available components H\ , H 2 as shown in Fig. 1.2. If the thickness of the layer d is 
much smaller than the cell size Ax then E t \^ and E t 2 can be both assumed to reside 
midway between H n and H t \ , therefore the ordinary FDTD update equations can be 
used. Then, the auxiliary variables X;y n , X tJy k are advanced in time and finally the 
new electric fields E n , Et 2 are obtained using the time domain equivalent of (1.32) 
and (1.33). 

Currently, the implementation of the above procedure in one dimension is in 
progress. 
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Chapter 2 


FDTD Predictions of Penetration 
Through Apertures in Conducting 
Enclosures 


The penetration of electromagnetic fields into conducting enclosures via apertures is 
an EMI issue that is relevant to all of aviation. The stories are numerous, of disrupted 
communications, disabled navigation equipment, and worse; due to the effects of EM 
sources external to the aircraft. 

Although there are many possible mechanisms for the penetration of fields into an 
aircraft (including direct penetration through composites, penetration through cracks 
and joints, conduction along cabling, etc.), it is usually the windows which admit the 
greatest. Consequently, it is of greatest importance to understand and be able to 
predict the field penetration through apertures. 

In this report, three cases of penetration are considered. The three cases consist 
of a modestly sized enclosure with three aperture configurations: one large aperture, 
two large apertures, and one small aperture. 

Contrary to initial intuitive estimations, these are not entirely trivial problems for 
the Finite-Difference, Time-Domain (FDTD) method. Particularly when the aperture 
is small relative to the size of the enclosure, the decay time of a pulsed excitation is 
enormous. However, there may be some advantages to introducing an artificial loss 
mechanism into the problem space to accelerate the decay of the late-time fields. 

Thank you to Dr. M. D. Deshpande of FDC/NYMA Langley Research Center, and 
Fred Beck and C. Cockrell of NASA Langley Research Center for providing the MoM 
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Case #1 



Figure 2.1: A 30 cm X 30 cm X 12 cm conducting box with a 20 cm X 3 cm aperture. 
A vertically polarized plane wave is normally incident on the aperture. 


predictions and measurements used for comparisons with the FDTD predictions. 


I. Case #1 

The enclosure is a rectangular box having the dimensions of 30 cm wide by 30 cm 
deep by 12 cm high. The first aperture case (Case #1) consists of a 20 cm wide by 
3 cm high aperture in the side of the box. The drawing in Figure 2.1 illustrates this 
small box with an aperture. 

As shown in the figure, the fields (a Gaussian pulse) are normally incident on the 
aperture. The Rayleigh number ( ( 3 ) for the pulse was chosen to yield a spectrum of 
up to 1 GHz. 

This geometry was discretized at a cell size of 2.5mm. This cell size is very 
small (A/120 at 1 GHz), but is necessary to resolve the aperture of Case #3. This 
discretization results in a solution space that measures 130 x 130 x 58 cells, for a total 
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Figure 2.2: The 30 X 30 X 12 cm box with the 20 X 3 cm aperture (Case #1) was 
discretized at 2.5 mm per cell (A/120 at 1 GHz). 

of 980,200 cells. A view of the mesh from the side of the box containing the aperture 
is pictured in Figure 2.2. Five pad cells separate the sides of the box from the six 
layers of PML. 

The pivotal issue in using a pulse excitation with the FDTD method is that the 
fields (or currents), particularly at the locations of interests, decay to zero. When 
the time-domain response has decayed to zero, it can be Fourier transformed into the 
frequency domain without error. In some cases, it may not be practical or possible 
to execute sufficient simulation time to achieve complete decay of the pulse response. 
In those cases, a reasonably good answer can be obtained by transforming the pulse 
response after it has decayed to a sufficiently small amount. However, the criterion 
for what is “sufficiently small” has not been formalized. 

After some “false starts,” Case #1 was executed for 131,000 time steps. This 
corresponds to a simulation time of 0.6 microseconds. In this 0.6/isec, the pulse 
could travel back and forth (across the long dimension) through the enclosure almost 
600 times. This simulation ran for approximately 46 (wall clock) hours on our SGI 
Octane. 

In Figure 2.3, the component of the time-domain electric field, at the center of the 
enclosure, that is copolarized with the excitation field is plotted with respect to the 
number of time steps. As shown in the inset, the field has decayed to an amplitude 
of 0.0002 V/m. As will be seen in Figure 2.4, this is sufficiently dose to zero to result 
in an accurate prediction for this case. 


(V/m) 





The results are presented here as elecric field shielding. The electric field shielding 
is a relation between the fields at the center of the box relative to the fields outside 
of the box. To realize the electric field shielding using the FDTD, the following 
procedure was performed. First, an empty FDTD solution space was created. This 
empty solution space does not need to be large. A space of 26 x 26 x 26 cells was 
chosen. The cell size of the empty solution space must match that of the problem 
solution space. Then the same pulse used in the problem simulation was launched in 
the empty solution space. The frequency response of the pulse at the center of the 
empty solution space was stored. Note that the number of time steps can be greatly 
reduced compared to that of the problem simulation. In this case, 32,000 time steps 
were used. This was more than necessary. Finally, the “electric field shielding” was 
computed according to 


electric field shielding = 201og( empty space frequency response 

problem space frequency response 

where “empty space frequency response” and “problem space frequency response” 
are the linear magnitude values of the z components of the frequency-domain electric 
fields computed at the centers of the empty and problem solution spaces, respectively. 
Note that with this definition, the fields at the center of the enclosure are highest 
when the electric field shielding is lowest (near 650 MHz for Case #1). 

The FDTD-predicted electric field shielding is plotted in Figure 2.4 with measured 
values and Method of Moments (MoM) predicted values. The MoM prediction model 
consists of the cavity mounted behind an infinite ground plane which includes the 
aperture. Thus the MoM prediction does not include any diffraction effects from the 
exterior edges of the enclosure. The correlation of the FDTD prediction is excellent 
compared to the MoM prediction and the somewhat noisy measurement. This noise 
is attributed to resonances in the screen room in which the measurements were per- 
formed [2]. 
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II. Case #2 


Case #2 consists of the same box as in Case #1, but with two 20 x 3 cm apertures. 
One aperture is at the same location as the one for Case #1, the second aperture is 
located in the center of the side of the box opposite the illumination. A drawing of 
the conducting enclosure with the two apertures of Case #2 is shown in Figure 2.5. 

All of the FDTD parameters (cell size, / 3 , number of pad cells, etc.) for Case #2 
are the same as those in Case #1. The predictions were again computed for 131,000 
time steps, and the same procedure was used to normalize the resulting frequency- 
domain field at the center of the enclosure relative to the input pulse response. 

The time-domain response of the pulse at the center of the cavity is shown in 
Figure 2.6. The field configuration and pulse response is different from that of Case 
#1 due to the different aperture configuration. The time-domain fields are somewhat 
lower for this reason. The peak response for Case #1 is 0.7262 V/m; for Case #2 
it is 0.4421 V/m. However, the amplitude of the late-time fields for Case #1 is 
approximately 1.93 x 10 -4 , or about 0.027% that of the peak. In Case #2, the 
amplitude of the late-time fields is about 4.8 x 10 -5 V/m, or about 0.011% of the 
peak. It is a reasonable conclusion that the fields inside the cavity of Case #2 have 
decayed more rapidly than those in Case #1 due to the additional loss mechanism of 
the second aperture. 

The only values available for comparison for this case are MoM predictions. One 
of the predictions is labeled “with internal coupling”, the other is labeled “without 
internal coupling” (with respect to the apertures). It is assumed by the author that 
these two predictions were intended to explore the effects of the coupling between 
the two apertures, and are not the most accurate MoM predictions possible. This 
is in consideration of the fact that this structure is slightly less demanding on the 
FDTD than that of Case #1, for which the agreement with MoM was excellent. The 
FDTD-predicted electric field shielding in Figure 2.7 is greater than that of the MoM 
predictions. However, the general shape is very similar, and the null occurs at ap- 
proximately the same frequency as in the MoM predictions. 
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III. Case #3 


Again, the enclosure for Case #3 has the same dimensions as those in Case and 
#2. This time, the single aperture is 0.5 cm in height and 10 cm wide. The vertically 
polarized pulse is again normally incident on the face of the box which contains the 
aperture, as shown in Figure 2.8. The geometry was again meshed at 2.5 mm per 
cell. The mesh of the aperture side of the box is shown in Figure 2.9. 

After executing Case #3 for 131,000 time steps, the time-domain electric field 
amplitude at the center of the box has only decayed to a value of 0.03 V/m, as seen 
in the inset of Figure 2.10. This high-Q cavity essentially represents a condition that 
is the opposite of Case #2 relative to Case #1: the smaller aperture of Case #3 
provides a smaller loss mechanism, and the fields decay very slowly. 

Is the late-time electric field amplitude of 0.03 V/m close enough to zero for the 
accurate transformation of the time-domain fields into the frequency domain? The 
electric field shielding was computed as before, and the results are plotted with MoM 
prediction and with measured values in Figure 2.11. The agreement between the 
measurement and the MoM prediction indicates that they are accurate, and that the 
FDTD prediction is significantly in error. The only agreement between the FDTD 
prediction, and the MoM prediction and measurement is the location of the null at 
700 MHz. Evidently, the time-domain fields have not decayed to a small enough value 
for the FFT to be accurate. 

One possible solution is to execute the simulation for a greater number of time 
steps. With the current configuration of the code, a radix-2 FFT is used to transform 
the time-domain fields into the frequency domain. That means that the computer 
memory reserved for the time-history of the fields to be transformed must be increased 
by powers of two: the smallest increase is double the current number of time steps. 
In this particular case, the array size could be increased from 131,072 to 232,177. 
Unless there is reason to believe that the late-time fields will decay to an acceptably 
small value within some intermediate number of time steps, the logical number of 
time steps to run is then 232,000. That would take at least twice as long to run as 
does the present number of time steps: approximately four days. Judging from the 
rate of decay evident in Figure 2.10, a doubling of the number of time steps is unlikely 
to be sufficient. This option is not particularly appealing. 
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Another possibility is to introduce some additional loss mechanism into the so- 
lution space. To enlarge the existing aperture or add an additional aperture would 
surely change the field distribution inside, thus invalidate the prediction. The walls of 
the enclosure could be made lossy; however, a third possibility is to artificially assign 
some small loss to the free-space cells of the solution space. This third approach is 
attempted. 

A. Case #3, artificial loss 

The FDTD prediction for Case #3 was computed again for 131,000 time steps, but 
with the free-space cells artificially assigned a conductivity of 0.000001 S/m. This 
change was implemented in the input (control) file. The second parameter of the first 
media property line defines the conductivity of the free-space cells. 

The time-domain electric field amplitude versus time step is shown in Figure 2.12. 
The late-time fields have decayed only slightly more than those for the lossless case. 
However, the electric field shielding has improved substantially as seen in Figure 2.13. 
The null at 700 MHz has increased to the MoM predicted value, and the erroneous 
peaks at approximately 410 MHz and 890 MHz have been reduced. However, the dis- 
crepancies between the FDTD prediction and the measurement and MoM prediction 
is still considerable, particularly below 300 MHz. 

Some improvement in the FDTD-predicted electric field shielding was observed 
in the previous case in which the free-space cells were assigned a very small value of 
conductivity. The obvious next step would be to increase this conductivity to obtain 
a more pronounced effect on the late-time fields. This case was repeated with the 
artificial conductivity of the free-space cells set to 0.0032 S/m. Furthermore, the 
number of time steps have been reduced to 16,000. 

In Figure 2.14, the time-domain fields at the center of the cavity are shown versus 
time steps, and the last 3,000 time steps are shown in the inset. The late-time fields 
have now been driven almost to zero within 16,00 time steps. 

A dramatic improvement in the comparison between the FDTD prediction and the 
measurement and MoM prediction is evident in Figure 2.15. However, the FDTD- 
predicted null at 700 MHz is slightly higher than both the measured and MoM- 
predicted null. 
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B. Case #3, variable artificial loss 

Finally, Case #3 is repeated after making a slight modification to the NEWS code 
to linearly increase the artificial conductivity assigned to the free-space cells. This 
variable conductivity is set to zero at time step zero, and increases to a value of 0.0064 
S/m at 16,000 time steps. The value of 0.0064 S/m was chosen to yield the same 
total loss experienced by the pulse: the integration of the instantaneous conductivity 
with respect to the simulation time. 

The time-domain field at the center of the cavity is plotted versus time steps in 
Figure 2.16. As seen in the inset, again the fields have been driven virtually to zero 
within 16,000 time steps. However, the difference in the rate of decay from that of 
the previous case is clearly evident. 

The electric field shielding for this case is compared with measurement and MoM 
prediction in Figure 2.17. The agreement is excellent. Note that the FDTD-predicted 
null at 700 MHz is now between that of the measurement and the MoM prediction. 

Why is the electric field shielding predicted by the FDTD slightly better when 
the artificial conductivity of the free-space cells is linearly increased than when it 
is set to a constant? It is speculated that there is an initial response to the pulse 
excitation that is critical to the accuracy of the final predicted spectrum. This initial 
pulse response is clearly seen in the time-domain field plots for Cases #1 and #2 
(within the first few thousand time steps), but is not obvious in Case #3. In Figures 
2.18 and 2.19 the time-domain field for the lossless Case #3 is compared with that 
using a constant artificial loss and a variable artificial loss, respectively. Notice that 
the time-domain response for the constant artificial loss case is different from that 
of the lossless case. However, the time-domain response of the linearly increasing 
artificial loss case is almost exactly the same as that of the lossless case for the first 
few thousand time steps. For the linearly increasing artificial loss case, the early-time 
loss is negligible and the initial pulse response is not significantly altered. The loss 
then increases to levels that are effective in causing a rapid decay of the late-time 
fields, enabling an accurate Fourier transformation into the frequency domain. 
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IV. Conclusions 


The FDTD-based NEWS code was used to compute the electric field shielding of a 
conducting box having three different aperture configurations. The FDTD predictions 
were compared with measurements and MoM predictions, where available. 

Due to the high Q of the cavities considered, the conventional second-order FDTD 
method required a very large number of time steps for a Gaussian pulse excitation to 
decay to a value that was sufficiently close to zero for accurate Fourier transformation 
of the time-domain fields into the frequency domain. With 131,000 time steps (in this 
case), however, the agreement of the FDTD predictions with measurements and MoM 
predictions were excellent. 

The loss mechanism which allows the time-domain fields to decay consists of the 
fields exiting the cavity through the aperture and terminating into the Absorbing 
Boundary Condition (ABC). When the aperture area is small relative to the size 
of the cavity, the late-time fields decay more slowly than when the aperture area is 
large. For small apertures, the number of time steps necessary to allow the late-time 
fields to decay can be enormous. An extremely large number of time steps has signif- 
icant drawbacks: long computation time, accumulation of dispersion, accumulation 
of numerical error, etc. 

This study suggests that the introduction of an artificial loss mechanism into the 
solution space may be a viable alternative to extremely large numbers of time steps. 
In this report, the free-space cells were assigned a low value of conductivity. This 
additional loss mechanism did force the transient fields to decay more rapidly. 

The speculations, based on the observations from this study, concerning the intro- 
duction of this artificial loss are as follows. When the free-space cells are assigned a 
very low value of conductivity, the loss experienced by the pulse is cumulative over a 
large number of time steps, and the greatest effect is on the late-time fields. When the 
free-space cells are assigned a higher value of conductivity, the initial pulse response 
of the cavity is altered, and the accuracy of the prediction is degraded. However, if 
the conductivity of the free-space cells is increased as a function of time steps from 
an initial value of zero, the initial pulse response of the cavity is virtually unaltered 
from its true free-space response. The artificial loss then increases to values that are 
effective in forcing the late-time fields to zero. 
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For Case #3 (a 30 x 30 x 12 cm box with a 10 x 0.5 cm aperture), 131,000 
time steps were found to be insufficient for the late-time fields to decay to suitable 
levels. By slightly altering the NEWS code to linearly increase the conductivity of 
the free-space cells from zero to 0.0064 S/m, the transient fields were driven to zero 
in 16,000 time steps. This resulted in a predicted electric field shielding that was 
in excellent agreement with measurement and MoM prediction, and was a reduction 
in computation time of approximately 87% over that of the unaltered NEWS code 
(which did not yet yield an accurate FDTD prediction). 

It should be noted that, strictly speaking, the use of time- varying material prop- 
erties violates an assumption made in the formulation of the conventional FDTD 
update equations. The update equations are derived using the time derivative of the 
electric field intensity, not the electric flux density. Therefore, there is an implicit 
assumption made that the permittivity (e) is constant with time. 

The introduction of a linearly increasing artificial loss to the free-space cells of an 
FDTD problem consisting of a cavity with an aperture was found to greatly reduce 
the computation time of the calculations, reduce the computer memory requirements, 
enable the solution of an otherwise possibly intractable problem, and to yield results 
that are in excellent agreement with other methods. However, additional cases should 
be examined to verify that this technique is generally applicable to this class of prob- 
lems. In addition, there may be advantages to other profiles (exponentially increasing, 
delayed, etc.) of the artificial loss that can be explored. 
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Case No. 1 



for Case between FDTD prediction, MoM prediction, and measurement. 
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Figure 2.6: The time-domain electric field at the center of the enclosure for Case #2. 
The fields have decayed to approximately .00005 V/m at the end of 131,000 time 
steDS. 




Case No. 2 



Figure 2.7: A comparison of the electric field shielding at the center of the enclosure 
for Case #2 between FDTD prediction, MoM prediction (with internal coupling be- 
tween the apertures), and MoM prediction (without internal coupling between the 
apertures). 
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Case #3 



E 



Figure 2.8: A 30 cm X 30 cm X 12 cm conducting box with one 10 cm X 0.5 cm 
aperture. 



Figure 2.9: The 30 X 30 X 12 cm box with the 10 X 0.5 cm aperture (Case #3) was 
discretized at 2.5 mm per cell (A/120 at 1 GHz). 
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Figure 2.10: The time-domain electric field at the center of the enclosure for Case 
#3. The fields have only decayed to 0.03 V/m at the end of 131,000 time steps. 






Case No. 3 



Figure 2.11: A comparison of the electric field shielding at the center of the enclosure 
for Case #3 between FDTD prediction (with no artificial loss mechanism), MoM 
prediction, and measurement. 


34 


(V/m 



0 26200 52400 78600 104800 131000 

Time Step 


Figure 2.12: The time-domain electric field at the center of the enclosure for Case 
#3 and with the free-space cells of the FDTD solution space artificially assigned a 
conductivity of 0.000001 S/m. The fields have decayed only slightly more than those 
of the lossless case after 131,000 time steps. 
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Case No. 3 



Figure 2.13: A comparison of the electric field shielding at the center of the enclo- 
sure for Case #3 between FDTD prediction (with a constant 0.000001 S/m artificial 
conductivity assigned to the free-space cells), MoM prediction, and measurement. 
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Figure 2.14: The time-domain electric field at the center of the enclosure for Case #3 
with the free-space cells of the FDTD solution space artificially assigned a conduc- 
tivity of 0.0032 S/m. The fields have decayed almost to zero after only 16,000 time 
steps. 
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Case No. 3 



Figure 2.15: A comparison of the electric field shielding at the center of the enclo- 
sure for Case #3 between FDTD prediction (with a constant 0.0032 S/m artificial 
conductivity assigned to the free-space cells), MoM prediction, and measurement. 
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Figure 2.16: The time-domain electric field at the center of the enclosure for Case 
#3 and with the free-space cells of the FDTD solution space artificially assigned a 
linearly increasing conductivity as a function of time steps. The fields have decayed 
to nearly zero after only 16,000 time steps. 
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Case No. 3 



Figure 2.17: A comparison of the electric field shielding at the center of the en- 
closure for Case #3 between FDTD prediction (with a linearly increasing artificial 
conductivity assigned to the free-space cells), MoM prediction, and measurement. 
The agreement is excellent. 
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Case 3 Time-Domain Fields 



Figure 2.18: A comparison of the initial time-domain electric fields at the center of 
the enclosure for Case #3. The solid line corresponds to the case of the free-space 
cells being assigned a constant artificial conductivity of 0.0032 S/m. The patterned 
line is for the lossless case. Note that differences between the two curves occur very 
early in the simulation. 
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Case 3 Time-Domain Fields 



Figure 2.19: A comparison of the initial time-domain electric fields at the center of 
the enclosure for Case #3. The solid line corresponds to the case of the free-space 
cells being assigned a linearly increasing artificial conductivity. The patterned line is 
for the lossless case. Note that there is little difference between the two curves for 
the first few thousand time steps. 
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Chapter 3 


Analysis of Cavity-Backed Slot 
Antennas: FDTD, FEM & 
Measurements 

I. Introduction 

Mutual coupling between cavity-backed slot (CBS) antennas mounted on a ground 
plane has been examined in the past using the Finite-Difference Time-Domain (FDTD) 
method and a hybrid Finite Element Method/Method of Moments (FEM/MoM) ap- 
proach. Predictions were compared with measurements performed in the Electro- 
Magnetic Anechoic Chamber (EMAC) facility at Arizona State University. In this 
report, simulation issues related to the modeling of CBS antennas with either FDTD 
or FEM/MoM are described. Furthermore, additional numerical results are presented 
and compared with measurements. 

Initially, FDTD modeling of CBS antennas is discussed, including an efficient 
source implementation for single elements and its extension to multiple elements. 
Then, the hybrid FEM/MoM approach is briefly described and its advantages are 
outlined. In the last section, FDTD and FEM/MoM are used to compute the input 
impedance of a CBS antenna and coupling between two identical CBS antennas. Fi- 
nally, some parametric studies of coupling are performed using the hybrid FEM/MoM 
code. 
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II. FDTD Modeling 

A. Source Implementation 

The computation of the input impedance of an antenna or the network parameters 
of a system of antennas involves the Fourier transform of the input voltages and 
currents. Therefore, using a transient excitation (pulse) the impedance or the network 
parameters can be determined over a frequency band by fast Fourier transforming 
(FFT) the time-domain data. The basic requirement for the FFT to work is to allow 
enough simulation time for the transient phenomena to decay. However, one of the 
main difficulties involved in FDTD simulations is that in some applications, e.g., 
resonant lossless structures, tens or even hundreds of time-steps may be required for 
the transient fields to decay. 

The voltage source that is used in FDTD simulations is commonly connected 
directly to the antenna; hard voltage source. This source is usually a pulse with 
a significantly greater than zero amplitude, only for a very short fraction of the 
total computation time. When the pulse amplitude drops essentially to zero, the 
source becomes effectively a short circuit. Consequently, any reflections from the 
antenna towards the source are totally reflected back. Additionally, the energy that 
is introduced in the computational space can be dissipated by either radiation, or 
absorption due to the presence of lossy media or lumped elements. For resonant 
structures, this radiation or absorption process may require a long time to dissipate 
the excitation energy. 

Many attempts have been recently made in order to reduce the required simulation 
time. A transient source for a microstrip line was introduced in [3]. This source 
is located near the FDTD outer boundary and when its amplitude falls to zero, 
the source is removed and replaced with a FDTD absorbing boundary. The main 
disadvantages of this method are that it can not be explicitly applied to arbitrary 
geometries, and the feed location should be far enough from any geometrical features 
so that no reflections return to the feed location before the source is replaced with 
the absorbing boundary. A similar approach has been applied to antennas that are 
fed by a coaxial cable. Instead of introducing an absorbing boundary during the 
calculation, a portion of the coaxial cable terminated in an absorbing boundary is 
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Source 




included in the FDTD calculation [4], [5]. It is reported in [5] that this approach 
is preferable to the gap source described in [6], since the fields decay more rapidly 
when the explicit coaxial cable is used. Furthermore, a more involved technique 
that reduces the simulation time is based on the use of signal processing methods to 
predict the voltages and currents at later times based on manipulation of computed 
results at earlier times. This technique [7], [8] proposes the truncation of the FDTD 
computations at a given time and use of the already calculated results to predict 
the ones at future times. However, the utilization of prediction methods increases 
the complexity of the FDTD calculation process (post-processing). Additionally, the 
prediction methods are complicated and their accuracy depends significantly on the 
choice of order of the prediction technique, which has to be made by the user and is 
not a trivial matter. 

A novel, effective and very simple technique to implement for reducing the simula- 
tion time is based on a source with an internal resistance that provides the excitation. 
Initially, this method was used in [9] to excite microstrip patch antennas. In addi- 
tion, the expression for a voltage source with an internal resistance in parallel with 
the free-space capacitance of the FDTD cell is given in [10]. However, the advantages 
of this method were illustrated and outlined explicitly only in [11]. 

Consider the excitation of an FDTD computational region with a voltage source 
at a certain mesh position (i s , j s , k a ). Figure 3.2 shows the equivalent circuit of a 
voltage source with an internal source resistance R s . The voltage at the output port 
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of the source, which is fed to the antenna, can be easily computed by applying Ohm’s 
Law to the circuit of Figure 3.2 


V mt = V,-I,R, (3.1) 

where I s is the current flowing through the source which is equal with the current at 
the input port of the antenna. This current can be computed using the integral form 
of Ampere’s law. Obviously, the input impedance of the antenna is calculated from 
the following equation: 



Figure 3.2: Voltage source with intrinsic resistance. 

The source resistance cannot be chosen very large, otherwise instabilities may 
occur due to neglecting the displacement current through the FDTD cell containing 
the source. This problem can be encountered by judicious choices of the values of the 
resistance or by taking into account the capacitance of the FDTD cell as discussed 
in [10]. A reasonable choice for R s is the value of the characteristic impedance of the 
transmission line; e.g., 50 ohms for a standard coaxial cable. 


B. Illustration of the effectiveness of the resistive voltage 
source 

In order to illustrate the effectiveness of the use of a voltage source with an internal 
resistance, it was decided to compute the input impedance of an air-filled cavity- 
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backed slot antenna analyzed in [12]. A three-dimensional (3-D) view of the cavity 
under consideration is demonstrated in Figure 3.3 and a detailed description of the 
geometry is shown in Figure 3.4. The input impedance of the cavity-backed slot 
antenna was measured in the Electromagnetic Anechoic Chamber facility at Arizona 
State University. In the experiment, the aperture antenna was mounted on a finite 
ground plane of dimensions 24 x24 in and the sharp edges were covered with absorbing 
material to reduce diffractions. Additionally, the antenna was rotated at an angle 
with respect to the principal axes and offset relative to the center of the ground, to 
direct the diffractions away from the aperture. Moreover, the input impedance of 
the same aperture mounted on an infinite ground plane was calculated in [12] using 
a hybridization of the Finite Element Method (FEM) with the Moment Method 
(MoM) and compared very well with measurements; showing that the dimensions of 
the ground plane do not have a profound effect on the value of the input impedance. 



Figure 3.3: A three-dimensional view of an air-filled rectangular cavity-backed slot 
antenna fed with a probe oriented in the y-direction. 

In the FDTD simulations, this aperture antenna was mounted on a 9 x 9 cm finite 
ground plane, which is smaller than the one used in the measurements to reduce 
the size of the computational space. In order to determine the input impedance of 
the antenna, the voltage and current at the feeding probe of the cavity-backed slot 
antenna have to be computed. The simulation time should allow the transient effects 
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Figure 3.4: A two-dimensional view of an air-filled rectangular cavity-backed slot fed 
with a probe oriented in the ^-direction. 


in both the voltage and current to decay so that they can be accurately Fourier 
transformed. As discussed previously, the voltage is a user defined function of time 
that decays quickly to zero, whereas the current may need a long time to decay. Two 
FDTD simulations were performed, and for both the cell size was 1.5 mm whereas the 
probe was excited with a voltage Rayleigh pulse. In the first case, the voltage source 
had zero internal resistance ( R s = 0), and in the second case, R s was set equal to 50 
ohms. Figs. 3.5 and 3.6 show the computed current at the probe exciting the cavity 
for the two cases, respectively. Obviously, the current in the case with no internal 
resistance ( R s = 0) did not decay to zero even after 64,000 time steps, indicating the 
resonant behavior and the high quality factor, Q , of the antenna. On the contrary, in 
the second simulation, the source current converged to zero very fast and the FDTD 
calculation time reduced significantly. 
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Figure 3.5: Current of the probe exciting the cavity-backed slot antenna for R s 
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Figure 3.6: Current of the probe exciting the cavity-backed slot antenna for R . 
50 ohms 




C. Calculation of the ABCD-parameters using a voltage 
source with internal resistance 


In the previous section a voltage source with an internal resistance was used to ef- 
ficiently compute the input impedance of an antenna. However, when a problem 
involves the calculation of the network parameters of a system of antennas, then a 
modified approach should be used. Here, because our analysis is related to mutual 
coupling between two antennas, only the evaluation of two-port network parameters 
will be considered. The computation time required for the calculation of the two- 
port network parameters that model a system of two antennas can be substantially 
reduced by using a voltage source with an internal resistance. The augmented system 
of two generic antennas along with the voltage sources is shown in Figure 3.7. This 
system can be thought of as the cascade connection of three two-port networks, as 
illustrated in Figure 3.7. 



Figure 3.7: Augmented system of two generic antennas with voltage sources. 

The proposed approach requires to initially compute the F-parameters of the 
entire system including the antennas and the load resistors. Then, the F-parameters 
of the system are converted to ABCD-parameters. It is known that the ABCD matrix 
of the cascade connection of two-port networks is equal to the product of the ABCD 
matrices representing the individual two ports. Note also, that the order of matrix 
multiplication must correspond to the order in which the networks are arranged, since 
matrix multiplication is not, in general, commutative. Therefore, the ABCD matrix 
of the entire system is equal to the product of the ABCD matrices of each sub-network 
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and can be written as 
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(3.3) 


where ABCDr s1 and ABCDr $2 are the ABCD matrices of the series resistances 
R s i and R S 2 , respectively, and ABCD ant is the ABCD matrix of the overall antenna 
system. The ABCD parameters of a. two-port network consisting of a series impedance 
Z between ports 1 and 2 (see Figure 3.8) are given by 

(3.4) 
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Figure 3.8: Two-port network consisting of a series impedance Z. 
Consequently, the ABCD matrices denoted as ABCDr s1 and ABC Dr s2 are given 


by 


Ar , j B Rs1 


1 R s i 

Cr s1 d Rs j J 


0 1 

Ar , 2 Br 3 2 


1 Rs2 

C R 32 Dr s2 _ 


0 1 
L 


(3.5) 


(3.6) 

whereas the ABCD matrix of the overall antenna system, ABCD ant , can be computed 
by the following expression: 
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(3.7) 


The computed ABCD an t matrix of the two antennas can be converted, if needed, to 
any of the other type of two-port network parameters using the appropriate conversion 
formulas. Following the method described, computation of the ABCD parameters of 
the two antennas can lead to great savings in the computational time. 
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III. FEM/MoM Hybrid Approach 


A vector finite element method (FEM) hybridized with the method of moments 
(MoM) is used in the analysis of cavity-backed slot (CBS) antennas. The FEM, 
which is based on linear tetrahedral elements, solves for the electric field distribu- 
tion inside the cavity or cavities, since we are interested in coupling. A spectral 
domain MoM is implemented through the continuity of the tangential magnetic field 
in the aperture(s) to solve for the field distribution in the exterior region; note that 
the cavity(ies) is(are) mounted on an infinite ground plane with a possible dielec- 
tric/magnetic overlay. The main drawback of the spectral domain MoM is that it 
becomes extremely slow with increasing the number of edges in the aperture(s). This 
problem is overcome by using an asymptotic extraction of the exponential behavior 
of the Green’s function; the asymptotic part is evaluated using a computationally 
efficient spatial domain integration. The excitation is based on either a plane wave 
incidence (scattering problem) or a coaxial feed model (radiation problem) imple- 
mented using the FEM. The radar cross section (RCS), input impedance, return loss, 
and gain versus frequency can be conveniently and efficiently calculated using the 
hybrid formulation. 

A. FEM/MoM formulation 

A two-dimensional (2-D) view of a CBS antenna mounted on an infinite ground plane 
coated with a dielectric layer is depicted in Figure 3.9. A vector FEM is implemented 
for the solution of the source-free Helmholtz’s equation 

V x ({fi,.}- 1 • V x E) - k 2 0 [e r ] E = 0 (3.8) 

inside the cavity volume (17) , where [e r ] and [/q] are, respectively, the permittivity and 
permeability tensors of the domain, and E is the unknown electric field. Full-tensor 
representation of t r and // r allows electromagnetic modeling of frequency dependent 
anisotropic materials. Dirichlet boundary conditions are imposed on all perfectly 
conducting surfaces, which implies that nxE = 0 on cavity walls. Using this boundary 
condition and the well-known Galerkin’s approach, the Helmholtz’s equation can be 
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written in a weak integral form given by 

JJJ^[tir)- 1 VxE)-(VxW)dV - kl JjJje r \-E-WdV 

+ jhvo //» • (W x a z ) clA = 0 (3.9) 

where W is the chosen testing function, k 0 is the free-space propagation constant, r/ 0 
is the free-space intrinsic impedance, and S denotes the area of the cavity aperture. 

For scattering, a linearly polarized plane wave, denoted by H mc , is incident on the 
aperture plane at an angle 0; with respect to the normal vector and <f>i with respect 
to the .T-axis. By imposing the boundary condition that the tangential magnetic field 
must be continuous across the aperture plane, one may write 

H = H !nc + H re/ + ^2 //_°° M(fc*, k y )- G (k x , k y , z\0) e j ^ x+k ^ dk x dk y (3.10) 

where H re -^ is the magnetic field reflected by the coated ground plane in the absence of 
the aperture, G is the spectral dyadic Green’s function for a coated conducting ground 
plane in the presence of a magnetic current source at the origin, and M is the Fourier 
transform of the magnetic current (E x a z ) just above the aperture. Substituting 
(3.10) into (3.9) and utilizing the inverse Fourier transform in conjunction with the 
definition that T = W x a z , the weak form of Helmholtz’s equation can be expresses 
as 

J /^(k] _1 V x E) ■ (V x W) dfl - kl J fjje r ] E • W dn 

- /£ M(k x , ky)- G (k X , ky, *10) • T (-k, -ky) dkjky 

= jkoVojjj, HiE + H^) • T dA (3.11) 

The finite element volume is discretized with tetrahedrons and the aperture with 
triangles. Thus, the electric field (E) inside the cavity is expanded in terms of a set 
of basis functions Wj’s and the magnetic current (M) in the aperture is expanded 
in terms of another set of basis functions T/s, where i corresponds to a global edge 
number. The second type of basis functions (TVs) were originally proposed by Rao et 
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al [13]. These are very similar to Wj’s with the only difference being the enforcement 
of normal continuity, rather than tangential continuity, across edges. The Fourier 
transform of these triangular basis functions is known in closed form [14]. 

The integral on the right-hand side of (3.11) represents the excitation vector which 
is evaluated using a pure spectral domain approach. The spectral integral on the left- 
hand side of (3.11) represents the admittance matrix for the exterior region of the 
cavity. The latter is evaluated using a mixed spectral/spatial domain MoM approach. 
The exponential behavior of the governing Green’s function is numerically extracted 
to improve the computational speed of the spectral integration. The asymptotic 
part is evaluated using a spatial integration which is known to be computationally 
more efficient. Thus, using the asymptotic extraction, the admittance matrix can be 
expressed in spectral domain as 

Yu = - 3 ~^ j fi(k x ,ky) • (g (k x ,k y ,z\ 0)- G h *|0)j 

r Tj( b y ) d JCgdk y 

- J ^JI^%(k x ,k y )-G h (k x ,k y ,z\0)-f j (-k x ,-k y )dk s dk y (3.12) 

where G^ is the dyadic Green’s function of a homogeneous space with e r = tf and 
Hr = (if ( e r an d f^r referred to the dielectric coating), and T; is the Fourier transform 
of the i th basis function. The first integral in (3.12), denoted as Vp , is evaluated using 
a pure spectral domain approach after converting to polar coordinates. The second 
integral in (3.12), denoted as , is evaluated using a spatial, instead of spectral, 
domain MoM 

V'f = -2k 2 f Ti(r )•[/ Tj(r') G h (r,r') dA] dA 

JA t UAj 

+ 2 [ V-Ti(r) [ Gh(r, r')V' • Tj(r') dA' dA (3.13) 

where A l (Aj) is the triangle supporting the i th ( j th ) basis function, r (r') is the 
position vector, and Gh is the Green’s function for a homogeneous medium with 
e r = ef and (i r = (if. 
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Once the unknowns (non-PEC edges) inside the cavity and aperture are assembled, 
a. linear system of equations is formed 


■ A/c/c M c/a 


' E c ' 
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M a / c M a ? a + V' a /“ 
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. ba . 


where \M) represents the finite element matrix given by 


Mij = |||(Vx W,) • H- 1 • (V x W,) cm - k] /// n W t • [e r ] • W,- cm (3.14) 


and [E] = [P] + \Y H ]\ the superscripts c and a, which denote cavity and aperture, 
respectively, are used to distinguish field interactions between the two regions. The 
rank of matrix [A/] is N c 4- N a where N c is the number of non-PEC edges inside the 
cavity, and N a is the number of non-PEC edges in the aperture. The rank of matrix 
[Y] is N a . 

The right-hand side vector, which is non-zero only in the aperture plane and 
expressed by the surface integral in (3.11), can be evaluated very conveniently using 
the spectral domain approach 
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for hard and soft polarization, respectively; i denotes the global number of an edge 
in the aperture, * indicates complex conjugate, and 


k xs = k 0 sin Qi cos <f>i 


k ys — A: 0 sin 9 t sin <$){. (3*17) 

The Green’s function definitions for Gf and G^ f as well as additional details in the 
derivation of (3.15) and (3.16) are provided in [12], [15]. 

For radiation, the formulation remains the same except for the excitation vector 
which is formulated in a different way. Instead of plane wave incidence, the antenna 
is now excited using a coaxial feed model which becomes part of the finite element 
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domain. The associated formulation for the coaxial feed model and the evaluation of 
the reflection coefficient at the coax/cavity interface are explicitly given in [12]. 

The far-zone radiated and/or scattered fields are calculated using the magnetic 
current distribution in the aperture. Specifically, these are given by 


Ee = 


jk 0 cos6e • ,fc ° r 


Gf(k xs ,k ys ) 


2tt 


Y Ej {-7 %{k xs , kys) sin 4 > + Tyjihs, kys) cos 4>} (3.18) 
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_ j k Q cos 0 e 

^ 4> — O Kjt 4> Si hys) 
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Y Ej k ys ) cos <f> + T*-{k xs , k ys ) sin </} (3.19) 

jeA 

where A is a triangle in the aperture. Knowing the far-zone fields, antenna charac- 
teristics such as RCS, directivity, gain and efficiency can be calculated. 


Z 



Figure 3.9: Two-dimensional view of a ferrite-loaded cavity-backed slot. 
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IV. Numerical Calculations 


Initially, in order to examine the accuracy of FDTD results computed using a source 
with an internal resistance, the input impedance of the cavity-backed slot antenna 
examined in Section II. was computed for three different cases. In all the cases, the 
feeding probe was excited by a voltage source with R s equal to 50 ohms. In the first 
case, the radius of the probe was not modeled (infinitely thin probe) and the cell 
size was 1.5 mm. In the second case, the radius was taken into account by using the 
thin-wire model and the cell size was 1.5 mm. In the third case, the cell size was 0.6 
mm and the probe itself was discretized along with the rest of the geometry. Fig. 3.10 
illustrates the computed input resistance and reactance of the aperture antenna for the 
three different cases. Also, the FDTD calculations are compared with measurements 
and with the results based on the hybrid FEM/MoM formulation which were reported 
in [12]. Obviously, the accuracy of the FDTD results depends greatly on the wire 
modeling of the probe that excites the antenna. Excellent agreement between the 
FDTD computations and the measurements is observed in the case where the probe 
was discretized. The improvement in accuracy for the last case can be attributed to 
the finer discretization and the enhanced modeling of the probe. 

Moreover, the coupling between two identical cavity-backed slot antennas (whose 
specifications are defined in Fig. 3.4) mounted on a square 9x9 cm ground plane, 
25 mm apart from each other, was computed by FDTD and FEM/MoM and the re- 
sults of the two methods compare very well (see Fig. 3.11 for the geometry specifica- 
tions and Fig. 3.12 for the coupling calculations). The discrepancies at the higher end 
of the band can be attributed to discretization errors. In these simulations the FEM 
mesh consisted of 75,874 elements with average edge size 0.16 cm whereas the FDTD 
used a cell size of 0.6 mm and a computational domain of 160 x 160 x 142 cells. Notice 
that two different FDTD methods were used. The first one, denoted as FDTD(2,2), 
represents a second-order accurate FDTD both in time and space. The second one, de- 
noted as FDTD(2,4), represents a second-order accurate in time and and fourth-order 
accurate in space FDTD. It can be observed that the FDTD(2,4) scheme does not 
provide a significant improvement in accuracy compared to the FDTD(2,2) scheme. 
However, it is expected that as the distance between the two apertures becomes large 
in terms of the wavelength, FDTD(2,4) will outperform FDTD(2,2). 
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Figure 3.10: Input impedance of a cavity-backed slot antenna. 
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Figure 3.11: Geometry of two identical cavity-backed slot antennas mounted on a 
square ground plane (for antenna specifications see Figure 1.4). 
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Figure 3.12: S-parameters of two identical cavity-backed slot antennas mounted on a 
square ground plane (for antenna specifications see Figure 1.4). 
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Parametric studies of coupling between two E-plane oriented cavity-backed slot an- 
tennas, shown in Fig. 3.13, were performed both numerically and experimentally; d de- 
notes the separation between the two antenna apertures. Only the hybrid FEM/MoM 
approach was used in these calculations because it was shown to be more efficient and 
accurate than FDTD. The advantages of FEM/MoM are attributed to the fact that 
this method does not discretize the space between the two cavities but rather the 
open-space interaction between the slots is accounted for in the MoM formulation. 
A sample of the results is illustrated in Figs. 3.14 and 3.15, where the computed 
coupling of the two antennas is plotted versus separation distance at different fre- 
quencies and compared with measurements. It is observed that the coupling between 
the two antennas diminishes at approximately the same rate as a function of aperture 
separation at all frequencies. Also, the agreement between the measurements and 
the FEM/MoM hybrid is excellent considering that the levels of coupling are below 
-20 dB. 



Figure 3.13: Top view of two identical cavity-backed slot antennas mounted on a 
ground plane (E-plane configuration). 

Finally, the FEM/MoM method was used to compute coupling versus separation 
distance for an echelon configuration of the two identical CBS antennas examined 
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Figure 3.14: Coupling versus distance for the E-plane configuration at: (a) 7 GHz 
and (b) 8 GHz. 








above (see Figure 3.16). The distance cl was fixed to 25 mm and the separation 5 
was varied. The calculations are illustrated in Figure 3.17 for different frequencies. 
For this type of configuration there exist no available measurements. Again, it is 
observed that the coupling drops as a function of distance in a similar way for different 
frequencies. Moreover, it can be seen that coupling drops with a considerably faster 
rate in the case of the echelon configuration than in the E-plane orientation. This 
coupling behavior agrees with what it has already been observed in patch antennas, 
and is due to the fact that along the ground plane the E-plane radiation pattern 
of an aperture antenna exhibits a maximum, whereas the IFplane radiation pattern 
exhibits a null. 



Figure 3.16: Top view of two identical cavity-backed slot antennas mounted on a 
ground plane (echelon configuration). 

V. Conclusions 

In this report an overview of the analysis of cavity-backed slot antennas in the context 
of FDTD and FEM/MoM was presented. Different numerical issues related to the 
modeling of such antennas were described. It was found that the use of a voltage 
source with an internal resistance in FDTD is indispensable for efficient computa- 
tions. Both FDTD and FEM/MoM were used to compute the input impedance of a 
single CBS antenna, and coupling between two elements, and compared very well with 
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Figure 3.17: Coupling versus distance for the echelon configuration. 

measurements. However, it was found that for coupling versus distance calculations, 
FEM/MoM is faster and more accurate than FDTD. This is attributed to the fact that 
FEM/MoM does not discretize the space between the two cavities but rather treats 
the open-space interaction within the MoM formulation. Therefore, it does not suffer 
from dispersion errors like FDTD which simulates the open space in addition to the 
cavity space. Also, the FEM/MoM computational space remains constant because 
only the CBS antennas need to be discretized in the hybrid methodology. On the 
contrary, FDTD has to also account for the space between the two antennas thereby 
yielding very large domains especially for large distances. Another advantage of the 
FEM/MoM approach is that it is more flexible in simulating arbitrary geometries 
without having the staircasing problem of FDTD. For example, coupling between ar- 
bitrarily oriented CBS antennas on a ground plane can be accurately computed using 
the existing FEM/MoM code, whereas the FDTD code would have to be modified to 
include contour-path methods in order to eliminate the staircasing errors. 
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Chapter 4 


Radiation Pattern Analysis: 
FDTD(2,2) versus FDTD(2,4) 

I. Introduction 

In our research lab, higher-order FDTD schemes were developed in the past and 
thoroughly investigated. Specifically, a second-order accurate in time and fourth- 
order in space FDTD scheme [FDTD(2,4)] was first examined. The characteristics 
of FDTD(2,4) were compared with the ones of the standard second-order in time 
and space FDTD [FDTD(2,2)]. It was shown that FDTD(2,4) exhibits significantly 
smaller dispersion than FDTD(2,2). Therefore, it is expected that FDTD(2,4) can 
handle electrically large domains more efficiently and accurately. However, there are 
some problems related to the behavior of FDTD(2,4) around discontinuities especially 
PEC ones. Some of these issues were examined in the previous report by analyzing 1-D 
domains and some ways to resolve them were presented. As far as truncation methods, 
the anisotropic PML technique was used to terminate FDTD(2,2) and FDTD(2,4) 
lattices. It was found that PML is a very effective truncation technique, and it 
performs equivalently on either FDTD(2,2) or FDTD(2,4) computational grids. 

After validating the FDTD(2,4) 3-D code along with PML, a more general FDTD (2,4) 
code was implemented. This code works in a similar way as the NEWS FDTD(2,2) 
code. Before the simulation starts, it reads a discretized geometry created by the 
Anastasia mesher and also an input file where all the problem parameters are speci- 
fied. This FDTD(2,4) code can handle any type of radiation problems like the NEWS 
FDTD(2,2) code. However, scattering simulation capability has not been imple- 
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merited yet. It should be emphasized that this is the first version of a 3-D FDTD(2,4) 
code, and therefore, this program is still under validation. This validation procedure 
consisted of different radiation problems, e.g., analysis of dipoles, monopoles, cavity- 
backed slot antennas, etc., where patterns, coupling and input impedances have been 
calculated. 

In this report, a few radiation problems are examined in order to investigate 
and illustrate possible advantages of FDTD(2,4) versus FDTD(2,2). Initially, the 
radiation patterns of a monopole antenna mounted on the tail of the NASA scale 
model helicopter are computed at 9.18 GHz. Moreover, the radiation patterns of 
a monopole antenna mounted on top of a rectangular box are calculated at 2 GHz. 
Both problems are analyzed using either FDTD(2,4) or FDTD(2,2) and the numerical 
results are compared with measurements. 

II. Numerical Analysis 

The first problem consists of a quarter- wavelength monopole at 9.18 GHz, mounted 
on the tail of the NASA scale-model helicopter. The geometry is shown in Figure 4.1. 
The principal radiation patterns of this monopole on the NASA helicopter were mea- 
sured in the Electromagnetic Anechoic Chamber facility at Arizona State University 
in the past. The physical length of the monopole was 8.17 mm (A/4 at 9.18 GHz), 
the cell size was 4 mm (A/8) and the computational domain was 184 x 440 x 98 cells 
or 21A x 51A x 1 1 A long. Both FDTD(2,2) and FDTD(2,4) were used to compute the 
principal plane patterns of the monopole. The numerical calculations are illustrated 
in Figs. 4. 2-4.4, where the predicted yaw, roll, and pitch plane patterns are compared 
with measurements. Considering these computations, it seems that FDTD(2,4) does 
not provide a significantly better accuracy than FDTD(2,2). The main reason for this 
is that the discretization errors of each method are mixed with staircasing errors. The 
artificial corners that have been created on the simulation geometry through stair- 
casing become diffraction points that affect substantially the shape of the radiation 
patterns. This effect becomes more dominant at high frequencies, as the one used in 
this simulation. Through a closer examination of the results, it seems though that 
FDTD(2,4) exhibits a slightly better accuracy than FDTD(2,2). This improved ac- 
curacy of FDTD(2,4) can be observed, for example, on the yaw pattern calculations, 
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and on the top, bottom and port side of the roll plane. 


Monopole 



Figure 4.1: Geometry of a monopole mounted on the tail of the NASA scale-model 
helicopter. 

Because of the coexistence of discretization and staircasing errors in the previous 
example, it was not clear which method is more accurate. Therefore, another problem, 
which is free of staircasing errors, was chosen to be analyzed. This problem consists 
of a monopole mounted on a rectangular box of dimensions 1 1 A x 5A x 2A at 2 GHz 
(see Figure 4.5). The physical length of the monopole was 60 mm (A/2.5), the cell 
size was chosen either 20 mm (A/7.5) or 5 mm (A/30), and the computational domain 
was 92 x 54 x 34 cells or 320 x 168 x 88, respectively. Both FDTD(2,2) and FDTD(2,4) 
were used to compute the principal plane patterns of the monopole. The FDTD(2,2) 
code was run for both cell sizes (20 mm and 5 mm) whereas the FDTD(2,4) code 
was run only for the coarser mesh. The predictions of the simulation with the finer 
discretization (5 mm or A/30) were used as a reference to determine the accuracy of all 
simulations performed with the coarser discretization (20 mm or A/7.5). The results 
of FDTD(2,2) and FDTD(2,4) for the coarser mesh are compared in Figs. 4.6 and 4.7 
with the FDTD(2,2) predictions obtained for the finer mesh. From these figures, it is 
obvious that even with such a coarse discretization of A/7.5, FDTD(2,4) still predicted 
very accurately all principal patterns, whereas FDTD(2,2) failed to predict some of the 
peaks and nulls in these patterns. This example clearly illustrates that FDTD(2,4) 
outperforms FDTD(2,2) in terms of accuracy. Furthermore, it is expected that as 
frequency increases, the accuracy of FDTD(2,2) will decline even further whereas 
FDTD(2,4) will still provide satisfactory results. Additional simulations at higher 
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Figure 4.2: Yaw plane radiation patterns of a monopole mounted on the tail of the 
NASA scale-model helicopter at 9.18 GHz. 
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Figure 4.3: Roll plane radiation patterns of a monopole mounted on the tail of the 

NASA scale-model helicopter at 9.18 GHz. 
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Figure 4.4: Pitch plane radiation patterns of a monopole mounted on the tail of the 
NASA scale-model helicopter at 9.18 GH^ 




frequencies were not performed as it was impossible to run a fine enough mesh such 
as A/30 at a frequency higher than 2 GHz with the available computational resources. 



760 mm 

Figure 4.5: Geometry of a monopole mounted on a rectangular box. 


III. Conclusions 

Radiation patterns of monopoles on helicopters and rectangular boxes were com- 
puted using FDTD(2,2) and FDTD(2,4). It was shown that FDTD(2,4) outperforms 
FDTD(2,2) in terms of accuracy and overall computational efficiency. 
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Figure 4.6: Roll plane radiation patterns of a monopole mounted on a rectangular 
box at 2 GHz. 
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Figure 4.7: Pitch plane radiation patterns of a monopole mounted on a rectangular 
box at 2 GHz. 
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Chapter 5 


A New Iterative Hybrid Method 

I. Introduction 

Analysis of slot antennas mounted on a flat conducting surface was performed in the 
past using a hybrid finite element /moment method (FE/MM) approach which is also 
referred to as the finite element/boundary integral (FE/BI) method. According to the 
FE/MM hybridization, the interior of the cavity is discretized using finite elements 
(tetrahedrons, prisms, bricks, etc.) whereas the exterior region is treated using a 
spectral or spatial domain moment method. As a result, the overall system matrix 
consists of the finite element matrix [A], which is highly sparse, and the admittance 
matrix [Y], which is dense; the two matrices are coupled through the unknowns in 
the aperture. 

For the analysis of a single cavity-backed slot antenna, it was shown in previous 
reports that the hybrid FE/MM approach has enormous advantages over the pure FE 
or FDTD methods. The reason for that is because both the FE and FDTD methods 
require discretization not only in the interior of the cavity but also in the exterior 
region. The discretization in the exterior region needs to be fine enough so that 
the absorbing boundary conditions (ABC’s), which are placed at some distance away 
from the slot, are accurate. This excessive discretization may become computational 
expensive, both in terms of CPU time and memory storage, since the number of 
unknowns could increase substantially. This computational load becomes even larger 
when coupling between two or multiple antennas needs to be computed. In such 
case, not only the interior of all slot antennas has to be discretized but also the free 
space surrounding them. Computer simulation of such a problem becomes almost 
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impossible using a pure FE or FDTD method. Even for the hybrid FE/MM approach, 
the computational burden becomes increasingly larger as the number of slot antennas 
increases. 

In this chapter, a new iterative hybrid FE/MM approach is proposed to compute 
efficiently and fast coupling among multiple cavity-backed slots that are placed ar- 
bitrarily with respect to each other and flushmounted on a fiat conducting surface. 
Electromagnetic coupling among these antennas is computed through an iteration 
algorithm. For simplicity, assume that there exist two cavity-backed slots spaced 
a given distance apart. For coupling computations, one antenna is excited with a 
constant voltage Vo whereas the other one is left open circuit. The new algorithm 
is as follows: at the first iteration, the antenna that is excited radiates in space in 
the absence of the second antenna. The radiated field induces a field in the aperture 
of the second antenna which propagates through its cavity. The second antenna, in 
turn, scatters the incident field in all directions thus affecting the field distribution 
in the aperture of the radiating slot. This perturbation in the field of the radiating 
element also changes the input impedance. The new field distribution is re-radiated 
toward the second antenna and the iteration process continues until convergence of 
the field in the two apertures is achieved. Once we reach convergence, the self and 
mutual impedances can be computed. 

The new algorithm becomes extremely advantageous when the two antenna con- 
figurations are identical; thus, reducing (as it will be discussed later in this chapter) 
memory requirements and CPU time. It also provides a mechanism to achieve a full 
hybridization of FE/MM with other numerical techniques such as the GTD, PO, or 
even FDTD. In this chapter, we will present the formulation of this iterative method 
and apply it to the problem of computing coupling parameters of two identical cavity- 
backed slots flushmounted on an infinite conducting ground plane. The results ob- 
tained using the proposed iterative method will be compared with results obtained 
using the direct FE/MM approach. A comparison of CPU times and memory storage 
will also be provided. 
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Antenna #2 


Figure 5.1: Geometry of two identical cavity-backed slot antennas mounted on an 
infinite ground plane. 

II. Formulation 

The problem under investigation is coupling among multiple cavity-backed slot anten- 
nas flushmounted on a conducting ground plane. For simplicity, we will first formulate 
the problem of coupling between two cavity-backed slot antennas similar to the ones 
shown in Figure 5.1. In order to compute the coupling parameters for such a configu- 
ration, one antenna is usually excited at the input terminal using a constant current 
7] whereas the input terminal of the second antenna is left as an open circuit. Ra- 
diation by the first antenna induces a voltage V 2 at the input terminal of the second 
antenna. Thus, the mutual impedance Z 21 can be calculated using 

Vo 

Z 21 = t (5.1) 

h 

.Similarly, the self impedance Zn is equal to 

v% 

Z u = (5.2) 

^1 
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Due to reciprocity, Z\i = Zi\ and Z 22 = Z\\. Knowing the complete set of Z- 
parameters, the corresponding S-parameters can be computed in a straightforward 
manner using standard transformation formulas found in microwave books. 

Using the hybrid FE/MM approach to solve this coupling problem, the interior 
of the two cavities, including the apertures, is discretized using tetrahedral elements. 
The free space in between the two slots will not be discretized but rather the inter- 
action between the two slots is taken into account through an integral formulation 
using the respective half-space Green’s function. This approach certainly provides an 
advantage over the FDTD and pure FE methods since they both require a brute force 
discretization in free space. By renumbering the unknowns inside the finite element 
volume so that the unknowns in the aperture appear first whereas the unknowns in- 
side the cavity appear last, the corresponding matrix system for a single radiating 
antenna is given by 


M a / a M a/c If E a 
A/ e / a M c/c E c 



(5.3) 


where the superscript a stands for aperture and c for cavity. The non-zero excitation 
vector b c is part of the finite element volume inside the cavity. In addition, the matrix 
M a/a is a pure method of moments matrix and is dense; M c t c is a pure finite element 
matrix and is highly sparse; the other two matrices M a ^ c and M c / a provide coupling 
between the field inside the cavity and the field in the aperture. Now, in case there 
are two cavity-backed slots in close proximity, the combined matrix system, after 
renumbering the unknowns, is given by 

M a/a M a,c j,a/a Q 1 [ £“ 1 0 ' 

M[ ,a M C J C 0 0 El b\ 

(5.4) 

7j/ B 0 Ml ,a M 2 q/c El 0 

0 0 Mt [ E 2 J 0 . 

The subscript 1 denotes antenna #1 and 2 denotes antenna # 2 . Note that for cou- 
pling calculations, when antenna #1 is excited, antenna #2 is left open circuit and 
vice versa. Thus, the right-hand side vector of the matrix system for antenna # 2 is 
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zero. The above matrix system can be written in a more compact and convenient 
form as follows: 


' Mi 

Tu ' 


'Ei - 


'hi ' 

. T 2l 

m 2 


. Ei . 


_ 0 _ 


Matrix M\ is the self matrix of antenna #1, M 2 is the self matrix of antenna #2, T 12 
and T 21 are the interaction matrices between slots 1 and 2, and b\ is the right-hand 
side vector corresponding to the excitation of antenna #1. It is important to mention 
here that the matrices Tu and T 21 are transpose of each other, thus only one of them 
is stored in memory. In case the two cavities/slots are identical, and therefor their 
mesh information must be identical, these two matrices are symmetric which means 
only half of the matrix entries need to be saved. 

As seen from the matrix system in (5.5), the electric fields E\ and E 2 are coupled 
only through the dense interaction matrices Tu and T 21 . If those interaction matrices 
■were to be zero, matrices Mi and M 2 would be completely decoupled, therefore, 
allowing us to solve for E\ and E 2 independently and more efficiently. However, since 
these matrices are not zero, we need to take into account the interaction between the 
two slots and the governing fields. 

The matrix system in (5.5) can be efficiently solved through an iteration algorithm. 
The iteration starts by setting the field E 2 to zero. This implies that antenna #1 
radiates in space as if antenna #2 was not present. In other words, for iteration k = 0 

E { 2 0) = 0 (5.6) 


and, therefore, we can write that 


Mi ■ E\ 0) + Tu • 4 0) = bi 

(5.7) 

Mi - E[ 0) = bi- Tu ■ 4 0) 

(5.8) 

> 

tq 

M ^ — s 

O 

11 

O- 

1 

O 

(5.9) 


where the superscript (0) indicates iteration k = 0. Thus, after the matrix system 
in (5.9) is solved, we can obtained E[°^ which simply represents the governing field 
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distribution of antenna # 1 in the absence of antenna #2. Then, from (5.5) we can 
write that 

M 2 • E ( 2 1] = 0 - T 21 ■ E\ 0) (5.10) 

where the product of T 2 \ • E^ requires N° x N% operations; jV“ is the number of 
unknowms in the aperture of antenna # 1 and N% is the number of unknowns in the 
aperture of antenna #2. Once the value of E 2 is updated, the iteration algorithm 
may continue on by updating E\ using 

M 1 ‘E[ 1) = b l -T n -E^ ) (5.11) 

The updating of the vectors E\ and E 2 will cont inue until an acceptable convergence 
tolerance is obtained. A measure of convergence can be defined based on the 2-norm 
residual given by 

/>. = Iff'-if-' 1 !, (5.12) 

Pi = ( 5 . 13 ) 

When both p\ and p 2 are less than a given tolerance, it means that convergence has 
been achieved. 

This hybrid method, which takes into account coupling between antennas through 
an iterative algorithm, is extremely advantageous when it is necessary to compute 
the coupling parameters of two identical antennas as a function of distance. For 
two identical cavity-backed slot antennas, only one geometry needs to be discretized 
since Mi — M 2 = M. Also, T 12 = T 21 = T which, in this case, is a square symmetric 
matrix. Thus, when coupling versus distance is calculated, the system matrix M does 
not change as a function of distance; the only matrix that changes versus distance is 
the interaction matrix T, which is relatively small (N a x N a ) and calculated using a 
spatial or spectral domain method of moments. As a result of this observation, the 
matrix M needs to be inverted only once, possibly using a sparse LU decomposition, 
whereas the updating of the fields E\ and E 2 is done relatively fast using a simple 
backward substitution. Even when the relative distance between the two antennas 
is changing, the matrix M does not need to be inverted again. Consequently, the 
coupling parameters as a function of distance may be computed with only minimum 


79 



computational effort which corresponds to the CPU time needed by the sparse LU 
algorithm to invert the matrix M. In addition, the memory requirements for storing 
the system matrix are reduced considerably, compared to a brute force discretization 
of the original two-antenna configuration. Memory savings also occur in storing other 
finite element arrays such as the connectivity information, nodal coordinates, edges, 
etc. 

More pronounced improvements in terms of CPU time and memory storage can be 
observed when the number of antennas is increased to three. In such a case, the self 
matrix M is identical for all radiating elements and, thus, it needs to be inverted only 
once for coupling versus distance calculations. In contrast, the number of interaction 
matrices increases from one to a total of three; i.e., Tn = T 21 , T 13 = T 31 , and I 23 = 
T 32 . Filling though these matrices as a function of distance is not computationally 
expensive, especially when using a spatial domain method of moments. 

III. Results 

In the previous section, the idea of field interaction between two distant radiating ele- 
ments through a novel iterative hybrid approach has been introduced and formulated 
for a generic problem. Each of the radiating objects is treated separately using the 
method of choice, whereas the coupling between the two objects is taken into account 
at subsequent iterations to correct the field distribution obtained from a previous 
iteration step. This approach is innovating and potentially powerful because it allows 
field interaction not only between two FEM domains, or an FEM and a. MM domain, 
but also between FEM and GTD or FEM and PO. 

Since this is the first time that this iterative approach is presented to the AHE 
program, it will be initially applied to the problem of coupling between two cavity- 
backed slot antennas mounted on an infinite ground plane, as is shown in Figure 5.2. 
The two cavity-backed slots are placed a distance D apart whereas the transverse 
dimensions of each slot correspond to those of an X-hand waveguide. The depth 
of the cavity is 7.7559 cm and the probe is located in the center of the horizontal 
dimension at a distance of 1.905 cm from the bottom face. The length of the probe 
is 0.6985 cm. 

The most obvious approach to solve for the coupling parameters of the two slots is 
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Figure 5.2: E-plane configuration of two identical cavity-backed slot antennas 

mounted on an infinite ground plane ( L = 2.286 cm, W = 1.016 cm). 

to perform a brute-force discretization in both cavities and apply the hybrid FE/MM 
approach. The interaction between the two slots is implicitly taken into account 
in the formulation. This approach has the disadvantage of creating twice as many 
unknowns, compared to the proposed iterative method, therefore requiring excessive 
memory storage and a long CPU time to solve the final matrix system. Besides, the 
matrix system needs to be re-solved every time the separation distance D between the 
slots is changed. In contrast, using the iterative approach, only one cavity needs to 
be discretized, thus saving memory storage and generating a smaller matrix system, 
which can be solved faster and more efficiently. Also, every time the separation 
distance is changed, the new field distribution can be conveniently computed using a 
simple back-substitution since the original self matrix is not a function of separation 
distance. This back-substitution can be performed provided the sparse LU or ILU 
factorization of the original matrix was obtained at the first iteration step. More 
substantial computational savings, both in terms of memory and CPU time, can 
be obtained when coupling among a larger number of identical cavity-backed slot 
antennas is computed. 

Using the iterative approach to calculate coupling, the two antennas are basically 
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treated separately whereas the field interaction between the slots is accounted for 
through an iteration procedure. It is not clear, however, how many iterations it 
takes for the field distribution in the two domains to actually converge. The first 
experiment performed in this study was to calculate the real and imaginary parts of 
the self and mutual impedances as a function of iteration number k at a frequency 
of 8 GHz and a separation distance of 1 cm along the E-plane. These results are 
illustrated in Figs. 5.3 and 5.4. For the sake of clarity, k — 0 indicates that the self 
impedance of antenna 1 was calculated assuming antenna 2 was not present. For 
k = 1, the field radiated by antenna 1 induced a field distribution on antenna 2 which 
was then scattered and altered the field distribution of antenna 1. This iteration 
procedure continues for a higher values of k. In both these figures, it is clear that the 
self and mutual impedances converge to a value at k = 3. For k = 0, which basically 
means that there is no interaction between the two antennas, the real part of the self 
impedance was found to be 37 ohms whereas the imaginary part was —86.15 ohms; 
the converged values, however, at k = 3 are 41.3 and —85.4 ohms, respectively. For a 
smaller separation distance D , these discrepancies between first and second iterations 
are even larger. Similar observations are made for the mutual impedance, too. For 
k = 0, the mutual impedance is absolute zero since the interaction has not yet been 
accounted for. As the iteration number increases, both real and imaginary parts of 
the mutual impedance converge to a constant value. 

This iterative technique was applied on the same antenna configuration, however 
with slot separation D = 2 cm, to calculate the self and mutual impedances as a 
function of frequency. The self impedance as a function of frequency is shown in 
Figure 5.5. The results obtained using the iterative technique for k = 5 are compared 
with the FE/MM direct approach which was validated against measurements and 
presented in Chapter 3. As illustrated, the comparison between the two techniques 
is excellent; note that the iterative technique is computationally more efficient. The 
mutual impedance as a function of frequency is shown in Figure 5.6. The agree- 
ment between the two methods is fairly good although there are some discrepancies 
for frequencies close to the parallel resonance of the slot. The reason is related to 
the accuracy of computing the coupling matrix between the two slots. This matrix 
is similar to a transfer function; knowing the input, the accuracy of the output is 
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Figure 5.3: Real and imaginary parts of the self impedance as a function of iteration 
number. The two slots are placed in an E-plane configuration with separation distance 
of 1 cm. The frequency of operation is 8 GHz. 



Figure 5.4: Real and imaginary parts of the mutual impedance as a function of iter- 
ation number. The two slots are placed in an E-plane configuration with separation 
distance of 1 cm. The frequency of operation is 8 GHz. 
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Mutual Resistance ( ohm ) Self Resistance (ohm ) 


solely determined by the accuracy of the transfer function. This coupling matrix was 
computed using a spatial domain method of moments which involves integration over 
triangular surfaces. These integrals are evaluated using a 13-point Gauss quadrature 
which is not strictly accurate, especially for large separations. In order to more ac- 
curately compute the mutual impedance, the 13-point Gauss quadrature needs to be 
improved by a higher-order integration. 




Figure 5.5: Real and imaginary parts of the self impedance as a function of frequency. 
The two slots are placed in an E-plane configuration with separation distance of 2 cm. 




Figure 5.6: Real and imaginary parts of the mutual impedance as a function of 
frequency. The two slots are placed in an E-plane configuration with separation 
distance of 2 cm. 

The return loss and mutual coupling as a function of frequency are depicted in 
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Figure 5.7: Return loss and mutual coupling as a function of frequency. The two slots 
are placed in an E-plane configuration with separation distance of 2 cm. 


Figure 5.7 for a separation distance of 2 cm along the E-plane. The iteration technique 
is compared with the hybrid FE/MM direct approach which is considered accurate 
enough. As illustrated, the agreement between the two methods is excellent. Even 
with one iteration the mutual coupling as a function of frequency compares very well. 
However, since the computational effort required for subsequent iterations is minimal, 
we often choose to compute the coupling parameters after 5 iterations have past. 

The most noticeable computational savings in terms of CPU time are realized 
when computing coupling parameters as a function of distance. The reason for achiev- 
ing such computational speed-up is because the system matrix is factorized (either 
using LU or ILU) only once. As the separation distance is varied, the updated field 
distribution inside the cavities is computed by using only a simple back-substitution 
which involves order of N operations, where N is the number of unknowns for the self 
matrix. Simply stated, the self matrix M is not a function of the distance between 
the two slots. 

The real and imaginary parts of the mutual impedance as a function of distance 
between the two slots are plotted in Figure 5.8. Three line traces are shown in 
each graph: one representing the direct approach and the other two representing 
the iterative approach. For the iterative approach, we chose to show the results for 
both “iteration 1” and “iteration 5” although there are no substantial discrepancies 
between the two traces. As those are compared with the direct approach, there exist 


85 




some small differences which are mostly attributed to the accuracy of the coupling 
matrix formulated through a spatial domain integral boundary equation. 




Figure 5.8: Real and imaginary parts of mutual impedance as a function of distance. 
The two slots are placed in an E-plane configuration and the frequency of operation 
is 8 GHz. 

Although Figure 5.8 illustrates some minor discrepancies between the direct and 
iterative hybrid techniques, as those were used to compute the real and imaginary 
parts of the mutual impedance, their comparison for mutual coupling (-S^) as a 
function of distance is very good. The corresponding graph is shown in Figure 5.9 
for a separation distance ranging between 30 and 160 mm. As illustrated, both 
methods give almost identical results. As far as the iterative technique is concerned, 
the difference in the predictions when using only 1 iteration versus when using 5 
iterations is more noticeable for relatively small separation distances. 

Before concluding this chapter, it is important that we discuss the computational 
requirements of the techniques used in this report to predict coupling of slot antennas. 
Basically, we have used two methods: a direct FE/MM hybridization and an iterative 
FE/MM hybridization. For the first method, the coupling mechanism between the two 
antennas is implicitly built into the matrix system, whereas for the second method, 
this coupling mechanism is explicitly taken into account through an iteration process. 
For the same discretization, the iterative method requires less than half the memory 
requirements compared to the direct method. Concerning computational time, it 
is good that we consider the case discussed in Figure 5.9. For this problem, the 
number of tetrahedrons used in the direct approach was 39,046 whereas the number 
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Figure 5.9: Mutual coupling as a function of distance. The two slots are placed in an 
E-plane configuration and the frequency of operation is 8 GHz. 






of unknowns was 41,120. For the iterative approach and using the same element 
mesh size, the number of tetrahedrons used in the computational domain was 19,543 
whereas the total number of unknowns was 20,583. Solving the matrix system for the 
FE/MM direct approach, a Conjugate Gradient Square (CGS) algorithm was used. 
The total CPU time spent to compute the mutual coupling between the two slots for a 
total of 14 points, as shown in Figure 5.9, was 288 minutes. For the FE/MM iterative 
approach, we had the choice of either using a sparse ILU decomposition to factor 
the matrix once and then use back-substitutions for the remaining iteration process, 
or using the CGS algorithm for all iterations. Both options were followed. When 
the sparse ILU decomposition was used, it took 134 minutes to obtain the matrix 
factorization; the remaining computations were performed within a very small time 
period (1-2 minutes, assuming 1 iteration per point). When the CGS was used, again 
assuming only 1 iteration per point, the total CPU time for 14 evaluation points was 
S3 minutes. However, as the number of iterations increases, the CPU time increases 
accordingly. For example, if 5 iterations are used, the total CPU time increases to 
265 minutes. 

As seen from these numbers, the FE/MM iterative approach results in both mem- 
ory savings and CPU speed-up. Specifically, if our main objective is to compute 
mutual coupling versus distance for a large number of points, the iterative approach 
using ILU factorization first and then back-substitutions might be the most attractive 
choice. On the other hand, if the number of evaluation points is small, the iterative 
approach using CGS might be preferred. 

IV. Conclusions 

In this chapter, we introduced a novel hybrid iterative approach to efficiently calcu- 
late coupling between two or more cavity-backed slot antennas mounted on an infinite 
ground plane. Each cavity is independently analyzed using the hybrid FE/MM ap- 
proach, whereas the field interaction between the antennas is taken into account 
through an iteration algorithm. This results in substantial computational savings 
both in terms of memory storage and CPU time. It was observed that this algorithm 
converges within only a few iterations. The results compare favorably with the hybrid 
FE/MM direct approach which has been used in the past to compute self and mutual 
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impedances of two cavity-backed slots. 

Another advantage of this new iterative approach is the potential for further hy- 
bridization of finite element with other methods such as the geometrical theory of 
diffraction (GTD) and/or the physical optics (PO). Field interaction between an an- 
tenna and an electrically large object in close proximity can be accounted for through 
this iterative approach in order to compute the overall radiation pattern as well as the 
input impedance. Most hybrid techniques only consider one-way interaction between 
the antenna and the object. In other words, the scattered fields by the object are 
not allowed to change the field distribution of the antenna and, therefore, its input 
impedance is the same as that one in the absence of the object. Using this innovative 
iterative algorithm, a full interaction between the two objects is possible. 
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Chapter 6 


Hybrid Techniques for Scattering 
and Radiation Problems 

I. Introduction 

Increasing the frequency of operation in helicopter communications provides a way of 
overcoming the limitations of line-of-sight communication. Satellite communication is 
an area from which modern helicopter technology is seeking for solutions to reach out 
over-the-horizon. In satellite communications usually higher than UHF frequencies 
are used, which results in the enormous size of computational domain. Therefore, it 
hinders use of a rigorous numerical method in the analysis of an antenna mounted on 
an airframe. As second choice, the adoption of hybrid techniques has been suggested. 

In fact, the FEM and PO were hybridized to analyze a monopole antenna on the 
APACHE helicopter in the UHF band. Its predictions for simple structures were 
provided along with measurements and exact solutions to validate the code. At the 
same time, the limitation of the method was also revealed; mainly due to having no 
account for diffractions and higher order interactions. Instead of using PO, the use of 
PTD allows inclusion of the diffractions in the analysis. However, the PTD solution 
is available only for simple structure; therefore, the implementation of PTD for a 
helicopter structure may take considerable amount of effort. Therefore, main focus of 
this report is on investigating and taking into account the higher order interactions 
through a hybrid technique. Here, higher order interaction means only the interaction 
between the antenna and its lit zone on the airframe. The diffractions can also be 
thought of as higher order interactions, but these are not included yet. 
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In this report, to get a better insight of the problem, a 2D hybrid code was devel- 
oped and used. As for the FEM analysis, more accurate boundary conditions than 
ABCs were implemented, referred to here as BIE (Boundary Integral Equation). Us- 
ing both truncation schemes, the FEM is hybridized with the PO which corresponds 
to the first order interactions. The reformulation was needed to account for higher 
than the first order interactions. In order to take the higher order interaction system- 
atically, an iterative algorithm is presented. It is also explained that the consideration 
for higher order interactions can be achieved using other than PO, e.g., GTD. The 
iterative algorithm shows that it allows the modification of equivalent currents at 
each iteration by having an updated right-hand side vector. 


II. 2D Radiation and Scattering 


The 2D-radiation problem can be decomposed into two decoupled TE Z and TM Z 
polarization analyses. For both cases, the scalar Helmholtz equation is given by 


d . dcj> s d . d(f> a 

~Tx la ai ) + ~a^ ia a^ ) + M - f 


where 
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for Ej-polarization and 
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( 6 . 1 ) 


( 6 . 2 ) 


(6.3) 


for /^.-polarization. The equation can be solved in conjunction with a boundary 
condition which can be derived from an ABCs or BIE [16]. For the finite element, 
triangular elements or quadrilateral elements are used to discretize the computational 
domain. If the triangular elements are used, a much larger number of unknowns will 
be generated compared to the quadrilateral elements. However, triangular elements 
do not require any computational effort to evaluate the integrals for the element 
matrix because all these integrals can be evaluated in an analytic way. In contrast, 
the quadrilateral elements lead to a smaller number of elements than the triangular 
elements, but the evaluation of the element matrix has to be performed in a numerical 
way. The truncation of the computational domain can be achieved by using ABCs. 
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The two ABCs implemented here are those of first (71, qi) and second order (72, q 2) 
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[VO'/p - M. 

This formulation is for the scattering problem; however, it can also be used for the 
radiation problem by setting < j> inc to zero. Based on this formulation, a 2D-FEM code 
was written. To validate the 2D-scattering code, scattering from a square cylinder 
was considered. The schematic of the problem is depicted in Figure 6.1. The FEM 
code does not have limitations in terms of coated materials on the scatterer, but a 
PEC square cylinder was chosen to compare with MOM solutions. Both triangular 
elements and quadrilateral finite elements were used in conjunction with the 2nd order 
ABCs. The mono-static echo width of a square cylinder predicted with the FEM is 
compared with the MOM results in Figure 6.2. The echo width can be computed 

from 

\ 6 SC \ 2 

= Jim 2^737 ( 6 - 9 ) 


Because of symmetry in the geometry, the echo width (<r) in other ranges is a repeating 
of Figure 6.2. As can be seen, the agreement between the FEM and MOM is very 
good, let alone the agreement between triangular and quadrilateral elements. The 
ABCs was applied at the distance of 0.35A away from the scatterer. It will also be 
interesting to test the accuracy of the 1st and 2nd order ABCs because the shape 
of ABC surface is almost square. Figure 6.3 shows that both the 1st and 2nd order 
ABCs work well for this problem, though all coefficients of the ABCs are derived from 
an assumption that the ABCs surface is a circular cylinder. 
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Figure 6.1: Scattering from a square cylinder when a plane wave is incident from an 
angle of <j) mc with respect to the x-axis. 

III. Boundary Integral Formulation 

When the FEM is used with ABCs to analyze an elongated geometry, there is always 
additional number of unknowns. However, the ABCs is a local boundary condition; 
therefore, the sparsity of the system matrix can be reserved. In other situation, such 
as two independent radiation systems are located close to each other, but it can- 
not be solved as one combined system, then ABCs may not be possible to apply 
because of the requirement of minimum distance. For this case, it is necessary to 
have a boundary condition that can be applied to the geometry as close as it can 
be. Having such a boundary condition, the overall radiation problem may be able 
to be represented through a coupling mechanism occurred between the two indepen- 
dent systems. Alternative boundary condition can be obtained from surface integral 
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Figure 6.2: Comparison of monostatic echo widths (|) calculated by the FEM and 
MOM. 

equations. Although, these are exact prior to discretization, they suffer from the fact 
that they are global conditions and create full matrices. However, if they are used 
just to obtain equivalent sources of the original source as an intermediate step, it may 
provide better conditions for the exterior region than ABCs. 

In the exterior region, the total field can be represented by the superposition of 
the original incident fields and scattered fields. Both are due to the source in the 
FEM domain, but it will be represented by the equivalent source. Therefore, the 
equivalent source may be arisen from an interior source or an interior scatterer. For 
the TM, polarization, these equivalent sources are expressed as 

1 f)F 

J z = z -fix ( lH t ) = H t = .. z - on dT (6.10) 

J KoV 071 

M z = i ■ ( hE t ) xii = E z on dT (6.11) 

For the TE Z polarization, the equivalent sources are given by 

J t = i • h x ( zH t ) = — H z on <9T (6-12) 
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Figure 6.3: Comparison of monostatic echo widths (|) calculated by FEM using 1st 
order ABC and the 2nd order ABCs. 


M z = z ■ (tH t ) x n = —E t = — ^ on dT (6.13) 

k 0 on 

where h is the outward unit normal vector to the surface, and i is the unit tangent 
vector defined so that h x i = z. In the exterior region, the equivalent source must 
satisfy EFIE for TM Z polarization or MFIE for TE Z polarization, which are given by 

E l z ne {p ) = M t (p) + z-V x [ M i {[))~H { Q ) (kR)dp + jh-j [ J z {p')^-H ( a\kR)dp' 

( 6 . 14 ) 

HT(p) = - i ■ V x [ J t (p')^H< 2 \kR)dp' + 3 — ( 

Jar 4 j n Jar 4 j 

(6.15) 

Let’s consider the TM Z polarization case shown in Figure 6.4. The interior region is 
discretized with a linear basis function of triangular elements. As consequence, the 
boundary T is broken into jVj short segments. The electric field on the boundary 
is piecewise-linear whereas the magnetic field is piecewise-constant; therefore, the 
basis function for EFIE requires the same type of basis function in order to provide 
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Figure 6.4: Hybrid between an FEM and an integral equation using linear basis 
functions on the boundary. 

continuity conditions for both fields. Moreover, using the simple testing function, such 
as Dirac delta functions located in the center of each cell, (6.14) can be expressed as 
a matrix equation 

e inc = Lk f + Mj* (6.16) 

where L and M are square matrix and e mc , k t and j z are all column vectors. The 
equation can be solved for j z to get 

j* = M-V nc - M _1 Lki (6.17) 

In fact, this is another type of radiation boundary condition and can be applied to 
FEM formulation instead of ABCs. Moreover, it is more general than the ABCs, 
since the shape of the boundary is not restricted. 
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In order to hybridized with the FEM, rewrite the discretized weak form of the 
Helmholtz equation in a matrix form 


In Ii6 o 


6j' n 


' 0 ' 

Ij,; E J 



— 

0 



. j* . 


L J 


(6.18) 


where the entries of 7 J)t -, 7,^ = 7^; and E represent interior interactions and can be 
obtained from 

[ [ -VNi-VNi-tferNiNidA (6.19) 

J Jr (l r 

and the entries of J can be expressed as 


Jij = -jhrj f BiBjdY (6.20) 

Ja r 

The basis function B, is the projection of Ni on the boundary, and B{ is the linear 
basis function for the electric equivalent current J z given by 

N b 

= ( 6;21 ) 


Finally, substituting (6.17) into (6.18) will lead to 

lii Iti> 

Im E-JM _1 L 

The combined system (6.22) forms a complete system for the solution of the total 
fields. Because of the global condition in the radiation condition, the part of the 
matrix with rows and columns associated with nodes on the boundary will be fully 
populated. The sparsity pattern of an example problem is depicted in Figure 6.5. 

It should be mentioned that care must be exercised in the computation of L and 
M matrix. Specifically, there are three cases, as can be seen in Figure 6.6, where 
the integrand becomes singular in every row, because the argument of 2D Green’s 
function becomes zero. To undertake more accurate evaluation, an asymptotic form 
of singular parts are subtracted from the original integral which results in an integral 
that can be computed by a numerical way. The remaining term which has to be 
added to the first integral contains the singularity but can be evaluated analytically. 
After calculating the total field, the equivalent current j z can be calculated from 
(6.17). Once the J and M are known, the scattered fields in the exterior region can 
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e 6 


— JM -1 e ,nc 


( 6 . 22 ) 
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be obtained using EFIE. In the far-held range, the scattered held can be calculated 
using the following formulae: 

E s = h+I 2 (6.23) 


h = - 


' Nb ■ f 

3 Z71 I 

n = 1 Jce 


Q jk 0 r’c° S (<p-v’) d Y' 


(6.24) 


l, ITT Nb r 

h = T\/4-E e ^/ B n co^-ip')t^ TC0 ^-^dT' (6.25) 

4 y 7T A.q xi— \ J cell 

The above formulation for TM^-polarization is implemented to predict the scat- 
tering width of an arbitrary shape of cylinder. To validate the BIE 2D code, the 
bistatic scattering width of a circular cylinder having 2A radius is predicted. Fig- 
ure 6.7 shows the bistatic echo width obtained by the BIE method in comparison 
with MOM results 
o 
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Figure 6.5: Matrix sparsity pattern after the hybridization. 
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Collocation (M) method with a linear basis function (N) 



I I Regions where the integrand becomes singular 


Figure 6.6: Cases where M and L elements become singular. 

IV. Hybridization of FEM with PO 

Although a hybridization of FEM with PO for 3D geometries was presented by the 
authors elsewhere, it was not a complete formulation. Three major issues pointed out 
to be improved are: 

a. Inclusion of higher order interactions between the source and the obstacle 

b. Augmentation of diffractions in the analysis 

c. Improvement of computation time needed for the PO integration 

With an aim to finding causes of those problems and hopefully provide solutions, it 
was necessary to simplify the problem in terms of size and complexity. Therefore, 
a good starting point was to develop an FEM/PO code in 2D-space. The geometry 
and coordinate system of the problem is shown in Figure 6.8. In the absence of any 
obstacle, the line source produces cylindrical waves. The exact solution for an electric 
line source in a cylindrical coordinate system is given by 

E z = -I^H^(hp) (6.26) 
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Plane wave Incident from <f>=0° 



Figure 6.7: Bistatic echo widths (^) predicted by BIE and MOM when the plane 
wave is incident from <p = 0°. 


H* = (6.27) 

Using the 2D FEM for TM Z polarization, the radiated field due to a line source in 
the absence of an obstacle can be calculated from (6.1) with / = — jk 0 Z o J z . By 
setting discretization size of A/20, the error of the fields for all observation angles was 
confined within 3%. In the presence of the obstacle, here a PEC body, the radiated 
fields will induce an electric current on the surface which may be evaluated using 
the PO approximation. Once the PO currents are available then, the scattered fields 
from the obstacle can be calculated from the EFIE with zero magnetic currents: 

E s = -jk oV [ Jpo(p')^-H { 0 2 \k 0 \p - p'\)dT' (6.28) 

Jht 4 J 

where Jpo is the PO current defined by 

J PO = 2(n x Hj ) if k ■ h < 0 (6.29) 
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The corresponding magnetic fields on the surface of the obstacle can be obtained 
using Maxwell equations: 

Ilf = - 'P-z x VE f (6.30) 

A/q 

Since the derivative of the 2D Green’s function w r ith respect to the observation points 
can be carried out with ease, the magnetic fields can be calculated without any ap- 
proximation. Owing to the PO approximation which is denoted as a forward reaction 
in the figure, the obstacle is replaced by the PO current in free-space; therefore, the 
exterior of the FEM domain must satisfy (6.14). At this time, E\ nc should be replaced 
by E s in (6.28). If the source were given, such as /(/?), the E s would be calculated 
exactly without the PO approximation using 


B'r(p) = - J / \p - A)dr‘ 


(6.31) 


It is also possible that other basis functions may be applied on the obstacle in order 
to get even better current approximation; however, they should be included in the 
matrix system, so their application will be limited only to small obstacles. Figure 6.9 
shows the backward reactions, which represent backward radiations from the PO 
source to the FEM domain and eventually perturb previous FEM solutions. The 
reaction between the forward and backward direction can be expressed as a matrix 
equation 


1-6 1 


I.'fc 

E - JM -1 L 
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u m 
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— JM -1 ( 


(6.32) 


where <5 m is a column vector having only one nonzero element at the node where 
the electric line source is located; whereas, e s is a column vector whose entries are 
scattered fields evaluated on the boundary nodes. Note that the second vector on the 
right-hand side account for the presence of an obstacle. In fact, the system matrix 
on the left-hand side and the first source vector on the right-hand side were already 
used to calculate the equivalent currents around the line source in the absence of the 
obstacle. Therefore, the matrix equation can be rewritten as 


Aj[ei + <5e] — Aie2 — bj -j- <5b 


(6.33) 


e 2 — ei + 6e = e! + A x *6b 


(6.34) 
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Note that the system matrix A will be the same as for the first step even for further 
iterations. From a computational efficiency point of view, LU factorization should 
be utilized to inverse the matrix A. The reason is that all that is required to obtain 
new 6e will be just backward and forward substitutions, if both L and U are stored 
from the first step. 


Lit region 


Obstacle 



Z 


* X 


Figure 6.8: Electric line source near an obstacle (forward reaction). 

The first example considered here is a line source near a strip. The distance 
between the source and the strip was 1.0 A and the width of the strip was 2.0 A. 
Figure 6.10 shows total fields computed by the FEM/PO in comparison to MOM 
predictions. The second order ABCs was used with the FEM/PO. The agreement 
between both results is shown to be excellent for both lit and shadow zones, even 
though the FEM/PO was performed without higher order iterations. From the width 
of the strip and the distance, the shadow region can be calculated in terms of ip angle 
as can be seen in the figure. It should be mentioned that the inclusion of higher oder 
effects will improve mainly the predictions in the lit zone. The improvement of the 
fields in the shadow zone can be accomplished only by supplementation of diffractions. 
Total field prediction for another example is shown in Figure 6.11. In this problem, 
the obstacle is chosen to be a circular cylinder having a radius of 1.0A. Based on 
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Radiation due to the current I (A/m) 

& 

Scattering due to the incident wave E b 



Figure 6.9: Electric line source near an obstacle (backward reaction). 

the geometry, the lit and shadow region can be defined as marked in Figure 6.11. 
The comparison to the exact solution exhibits good agreement in most of the lit zone 
except near boundary called transition regions. Again, the higher order interactions 
are not yet included, since those are not incorporated into our 2D FEM/PO code. 
Instead of using ABCs to calculate the equivalent source, the 2D-BIE code is used 
to calculate PO currents. The comparison between the BIE/PO and ABC/PO is 
depicted in Figure 6.12. Basically, the total fields predicted by both scheme were 
identical. However, the BIE/PO will be used as for the future analysis, because it 
saves lots of computation efforts when higher oder interactions are taking into account. 
The main reason is its use of collocation method. Because of that, the calculations 
of the interactions are needed only at the matching points. 

V. Hybridization of FEM with GTD or PTD 

For the case where the obstacle has smooth and flat surface, the FEM/PO works well. 
If the body has more lit region from the source point of view, the performance of the 
FEM/PO would be even better. However, for the special geometry such as wedges, 
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Figure 6.10: Total radiation fields due to a line source near a strip. The distance 
between the source and plate is 1.0A and the width of the plate is 2.0A. Comparison 
between FEM/PO and MOM predictions 


the FEM/PO may fail but GTD solutions may be valid. The same approach used for 
the hybridization between the FEM and PO can be expanded for GTD as well. The 
problem under consideration is shown in Figures 6.13 and 6.14. In fact, the exact 
modal solutions are available in the form of an infinite series [17]. 


F t _ pi , ™ _ / a v J„(pp)H ( J\pp') sin[i>(^ - a)] sm[v(<f> - a)] p < p‘ 
z z 2 \ E„ Gt ’JviPp^H^KPp) sin[v(^ - a)] sinj v{<j> -a)) p> p' 

(6.35) 

where a, and v are 


_ —Qjpirle 
2(?r — a) 


(6.36) 


mu 

2(7 r — a) 


(6.37) 


The magnetic field components can be obtained by using the Maxwell’s equation as 


rt _ - 1 j £v va v J l/ {0p)H^ ) {^p') sin[u(<£' - a)] cos [v{<f> - a)] p < p' .. 

p jtjjpp \ Hv va vJv{Pp')H[ 2 \l3p) sin[v(^' - a)]cos[u(^ - a)] p > p' 
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Figure 6.11: Total radiation fields due to a line source near a cylinder. The distance 
between the source and the cylinder is 2.0A and the radius of the cylinder is 1.0A. 
Comparison between FEM/PO predictions and the exact solution. 


= A. / £» sinbW' - a)] sinM* - a)] p< p' 

* 1 sin[i;(0' - a)] -siri - rt)] p>p’ 

Using the above total magnetic fields, the surface currents on the wedge can be 
obtained by 

Jgtd - n x (pH], + 4>Hl) = (n x H l p ) z (6.40) 


Instead of using Jpo in (6.28), Jgtd can be substituted to get scattered fields which 
will allow us to calculate higher order interactions in the FEM analysis. However, 
if both source and observation points are in the far-held range at the frequency of 
operation, then the asymptotic solution or GTD should be used because of its high 
convergence rate and accuracy. In GTD, the total scattered fields consist of GO fields 
and diffracted fields. At the given source point and observation points, GO fields can 
be calculated easily using geometrical optics; whereas diffractions fields are obtained 
using the diffraction coefficients which are available in [17]. 

The physical theory of diffraction (PTD) for the wedge can also be applied in 
conjunction with the FEM. Based on the PTD theory [18], the PO current needs to 
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Figure 6.12: Total radiation fields due to a line source near a cylinder. The dis- 
tance between the source and the cylinder is |A and the radius of the cylinder is |A. 
Comparison among FEM/PO, BIE/PO predictions and the exact solution. 

be supplemented by the contributions from the nonuniform current 

Km = Epo + (El - E d P0 ) = E™' + £"* (6.41) 

where the E represents diffracted fields calculated from Keller’s diffraction coeffi- 
cients. The electric field due to the physical optics current on the illuminated face 
is already calculated, so here only the diffracted field due to the PO current will be 
addressed. Based on the method of stationary phase, the fields due to the PO current 
is referred to as the end-point contribution. Using the geometry shown in Figure 6.13, 
the nonuniform current [19] is known to be 

• * — jx^ 

Jnu = -2Jpo l F{x)- e —~ (6.42) 

7T [ Zjx 

where F(z) is the Fresnel function and defined by 

poo . 3 

F(x) = / e~ jt dt (6.43) 

J X 

More general expressions for other than the half-plane are also available. In order to 
treat the obstacle using the PTD, the approximated PO current in (6.28) should be 
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replaced by the PTD currents given by 

Jptd = Jpo + Jnu (6.44) 

If the obstacle is a cylindrical body, the UTD can provide the scattered field not 
only in the lit and shadow regions but also in the transition regions without any 
discontinuity. Figure 6.14 shows a detailed geometry for scattering from a cylinder. 
For an observation point in the lit region, the reflected field is given by 

E'(p) = E*(Q')R, (6.45) 

where E‘(Q T ) is the incident field at the reflection point at Q r , and here only TM Z or 
soft-polarization is considered. Distance parameters are defined in the Figure 6.14. 
The generalized reflection coefficient for this polarization is given by 

*■ = (S 11 - F(Xr)] + M (6 - 46) 

where 


a. F(x ) is the transition function already used in connection with edge diffraction 
in GTD and PTD. 

b. The argument of the transition function is X p = 2 k 0 L p cos 2 (0;) and the distance 
parameter L p is given by 

~T A 

(6.47) 


^ s r + s' 

c. The (p is the Fock parameter associated with the reflected field in the lit region 


C p = — 2m(Q r ) cos 6 t 


(6.48) 


d. The value of the curvature parameter at the reflection point is expressed as a 
function of the radius of curvature a 0 of the surface at Q r 

m(Q,)= (6.49) 
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e. Finally, the function P s is the Pekeris caret function and can be written in terms 
of Fock scattering function p*(x) [20]. 

, . , e -3*P 

P,(x) = p(x)e (6.50) 

In the shadow region, the diffracted field is given by 

E d M = E\Q’)T, e -Ff- (6.51) 

E l (Q f ) is a GO type incident field at the attachment, point Q\ Observation point pd 
may be in the shadow region and even on the shadow boundary. The UTD surface 
diffraction coefficient, calculates the surface waves launched from two detachment 
points on the circular cylinder. Figure 6.15 shows decomposed GO field and the 
diffracted fields due to a line source in the presence of a circular cylinder. The 
distance between the source and the cylinder was 2.0A and the radius of the cylinder 
was 1.0A. As can be seen, the GO fields are zero in the shadow zone, so that the 
diffractions will be the main contribution for the interactions. 

VI. Conclusions 

In this report, an iterative algorithm was developed to analyze the radiation due to 
a source in the presence of a large body. In order to take the presence of an obstacle 
into account, the right-hand side vector of the system equation is modified at each 
iteration. Updated right-hand side vector contains scattering fields which represent 
backward radiation due to an induced current on an obstacle. To get a better insight 
of the problem, a 2D hybrid code was developed and used. As for the FEM analysis, 
BIE formulation was added to the conventional mesh truncation scheme, ABCs. 
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Figure 6.13: FEM analysis in the presence of wedge. 
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Figure 6.14: FEM analysis in the presence of a circular cylinder. 
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Figure 6.15: Decomposed GO fields and diffracted fields due to a line source in the 
presence of a circular cylinder. The distance between the source and the cylinder was 
2.0A and the radius of the cylinder w r as 1.0A. 
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Chapter 7 

Spectral Methods 

I. Introduction 

In the past, the authors introduced a family of spectral methods to provide an al- 
ternative way to handle efficiently electrically large computational domains, such 
as helicopter airframes. Spectral methods are very powerful due to their high accu- 
racy, exponential convergence, and negligible dispersion and dissipation. All these at- 
tributes are extremely important in computational electromagnetic problems. Some 
of the key elements in the formulation of spectral methods were presented by the 
authors elsewhere. The high accuracy of the spectral methods was demonstrated 
through numerical experiments and compared to the accuracy of the finite-difference 
methods. 

In this report, the spectral methods are applied for the first time to the solution 
of one-dimensional electromagnetic problems. Here, the ability of spectral methods 
to solve electromagnetic problems will be demonstrated, and their accuracy will be 
illustrated through numerical experiments. Their accuracy will also be compared to 
the one of the standard FDTD algorithm. Spectral methods, as expected, will exhibit 
orders of magnitude better accuracy than FDTD because of their negligible disper- 
sion. Different types of numerical experiments are presented here. Initially, the entire 
domain Chebychev collocation method is applied to solve an 1-D TAP-mode propa- 
gation problem. Furthermore, the PML technique is used in the context of spectral 
methods to simulate “open space” problems. Then, the ideas of interior domain grids 
and fictitious points are presented and applied to a multi-domain approach in the last 
two sections. This material is based on a recently published paper [21]. 


Ill 


From the material presented here, it will be shown that the spectral methods are 
very promising for numerical electromagnetic analysis. Our ultimate goal is to apply 
and extend these methods initially to generic 2-D and 3-D problems and then to 
radiation problems. 


II. Solving Maxwell’s Equations 


Here, the one-dimensional (1-D) Maxwell’s equations for the TAP mode given by 
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(7.1) 

- <?e}\ 

(7.2) 


are examined. Direct application of the entire domain Chebychev collocation method 
is quite straightforward. The spatial derivatives are approximated by using the dif- 
ferentiation matrices (DM) and an appropriate grid. Assuming a Chebychev grid 
Xi, i = 0, 1, ..., N and assuming grid values for the fields E Zi and H Vt , i — 0, 1, ..., N 
we can write the semi-discrete form of Maxwell’s equations as: 


Pi IT 1 

-w = j ( 0£ > - *<) < 7 - 3 > 

^ = ‘ (DR, - oE.) (7.4) 

where H_ y and E* a r^ the vectors containing the grid values of the field components 
and D the differentiation matrix. 

The first numerical experiment simulates a 1-D domain, one wavelength long, for 
the T M z - mode. The domain was excited with a sinusoidal source at the left hand side, 
and was terminated at the right boundary using the exact solution. This problem 
was solved using either the Chebychev collocation method or the FDTD second- 
order accurate scheme for different number of points. The error was computed in the 
maximum norm at time t = 2T, where T is one period. Figure 7.1 compares the 
FDTD results with the ones of the spectral method. Obviously, the spectral method 
outperforms the FDTD method in terms of accuracy. It is quite impressive that in 
order to get the same accuracy as the one that the spectral method provides with 
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a computational domain of 70 points, an FDTD domain of 1024 points is needed. 
Therefore, by using spectral methods, high accuracy can be achieved without very 
fine spatial discretization. The main drawback though of spectral methods is that 
they exhibit very restrictive stability criteria, forcing us to discretize very finely in 
time. However, there have been proposed approaches that deal successfully with this 
problem and accomplish to relax the time-step restrictions. 



Figure 7.1: Maximum error in the solution of the 1-D TM mode at t=2T. 


One main concern in electromagnetic simulations is the solution of unbounded 
space problems, often called “open-space” problems. The correct treatment of such 
problems is to simulate a computational domain large enough to enclose the structure 
of interest, and truncate it by applying appropriate artificial boundary conditions on 
the outer perimeter of the domain, which simulate its extension to infinity. These 
boundary conditions should allow the outward propagating wave to exit the domain 
and suppress spurious reflections of this wave to an acceptable level. Depending 
upon their theoretical basis, outer grid boundary conditions of this type have been 


113 




called either radiation boundary conditions (RBCs) or absorbing boundary condi- 
tions (ABCs). Initially, ABCs were obtained by applying the theory of one-way wave 
equations derived by Engquist and Majda [22]. Recently, a novel and pioneering ABC 
was introduced by Berenger [23], who used a non-physical absorbing (lossy) mate- 
rial adjacent to the computational boundary to truncate the domain, and he called 
this approach Perfectly Matched Layer (PML). The PML exhibits characteristics that 
permit electromagnetic waves of arbitrary frequency and angle of incidence to be ab- 
sorbed while maintaining the impedance and velocity of a lossless dielectric. Berenger 
reported reflection coefficients for PML in two dimensions significantly better than 
second- and third-order one-way wave equation (OWWE) based ABCs. 

Here, the PML concept is applied in the context of spectral methods to truncate 
an 1-D domain. According to Berenger the PML material should have properties so 
that the following relation is satisfied: 

a a* 
e /.i 

where e and ft are the properties of the material we want to match, e.g., in the case 
of free space e = e 0 and // = // 0 . To implement the PML we just need to assign to 
some of the grid points towards the end of the domain such material properties. Two 
are the key elements in an effective PML implementation: 

• Define a desirable and acceptable reflection coefficient. 

• Define an appropriate conductivity grading so that no rapid changes in conduc- 
tivities occur in the domain creating numerical reflections. 

The effectiveness of PML is shown here by solving again an 1-D domain which is 
terminated by PML with a desirable reflection coefficient 10 -6 and backed by PEC 
walls. The grading of the conductivity is quadratic. The Chebychev collocation 
method again is used to solve a half-wavelength long domain. The right hand side 
end of the computational space is the PML, and a Gaussian pulse was excited at 
the left hand side end. The pulse propagates along the domain (towards the right) 
until it enters the PML where is gets absorbed. Many simulations were performed 
with different number of grid points for the PML, and it was found that PML works 
equivalently well in the spectral methods as in the FD methods. A sample of these 
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results is shown in Figure 7.2 where frames of the propagating pulse are plotted in the 
same graph for a case where 16 points were used for the PML termination. Clearly, 
this figure illustrates the absorption of the pulse in the PML. Therefore, PML can 
be used to terminate effectively computational domains in the context of spectral 
methods. 



xlX 

Figure 7.2: Traveling pulse in 1-D domain terminated with PML. 


III. Interior Domain Grids and Fictitious Points 

The material that is illustrated here is based on a recently published paper [21]. As 
claimed in this paper, “Accuracy can be gained by clustering a grid less and using the 
approximations only well within the interior”. Given a grid of points j = 1, ..., N, 
the interpolation error for a smooth function is given by: 

f( N )(£) N fW(£) 

/(x) - pn-i(x) = ^ no* - *;) = $ ( x ) ( ? - 6 ) 
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for some £ £ [xi,xjv]. The ip(x) is called the remainder, and it is known from 
interpolation theory that it is minimized in the max-norm over [-1,1] when the nodes 
are forming a Chebychev grid. Figure 7.3 illustrates the remainder for a Chebychev 
and an equispaced grid for N = 10. Notice that both grids are interior meaning 
that they do not contain the ends of the interval. Notice that the remainder of the 
equispaced grid is large near the ends of the interval (Runge phenomenon) whereas 
the remainder of the Chebychev is larger over most of the interval. Thus, maybe a 
better choice for a grid is to have node distribution between the Chebychev and the 
equispaced ones. 

Other node distributions are defined in [21]. According to this paper a closed grid, 
meaning that the endpoints are included, is defined implicitly by 

jv = £T (i -* J r , «fc (7.7) 

c _ r (§-0 
1 v^FT (1-7) 

where Xj for j = 0,1,..., A r are the grid points. However, an open grid, i.e., not 
containing the end points, is defined by 

= £' c,(i - i ! r<ix (7.9) 

where Xj for j = 1 ,...,jV are the grid points. These equations for either the open 
or the closed grids can be solved numerically for a specific 7 < 1 to obtain the grid 
points. Notice that 7 = 0 corresponds to the equispaced grid and 7 = 0.5 corresponds 
to the Chebychev grid. Even though the optimal choice of grid in terms of error is 
the Chebychev one (7 = 0.5), it was suggested in [21] that a different choice of 
0 < 7 < 0.5 can possibly have some advantages. To show how the locations of the 
nodes are changing, depending on the choice of 7, the equations corresponding to 
the open grid were solved for N = 20 and for different 7, The results are plotted in 
Figure 7.4 where each contour illustrates the positions of each node versus 7. 

Another idea used in [21] is the concept of fictitious points, which is simply to 
augment a grid with additional nodes at which the function values are not directly 
known nor updated. This yields an interpolant that is used only in the interior of the 
extended interval. To illustrate this, we construct a grid by taking a 10-node open grid 
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for 7 = 0.35 and by adding to it one fictitious point per side of the interval by reflecting 
the outer points about ±1. Then, we plot the interpolation remainder for this grid 
and compare it with the one of the corresponding Chebychev grid (N = 10 and open 
grid) in Figure 7.5. It is seen that the fictitious point grid has a decreased error even 
though the minimum spacing between true points has increased by approximately a 
factor of 1.5. 



Figure 7.3: Interpolation remainder for N = 10. 


IV. Multi-domain 

In this section the multi-domain approach proposed by Driscoll in [21] is presented. 
The formulation will not be presented here but only briefly discussed. Basically, 
this multi-domain method breaks the domain under analysis into a number of sub- 
domains. The grid of each sub-domain is constructed by an interior node distribution 
(open grid) and m fictitious points determined by reflection of grid points about 
the boundaries of the sub-domain. The values of the fields at the fictitious points 
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are determined by enforcing certain conditions, such as the continuity of the fields, 
their first time derivative, etc., or by enforcing the boundary conditions known by 
electromagnetic theory in case the boundary occurs between two different media. 
Details of the approach are described in [21]. 

The multi-domain approach has the advantage of using smaller differentiation 
matrices and a larger time-step because each sub-domain is formed by a number of 
points considerably smaller than the one of the entire domain method. Note that in 
all the numerical simulations illustrated in this section one fictitious point per side of 
a subdomain is used. 

The first numerical experiment was based on a three sub-domain code which was 
written for validation. The entire free-space computational domain was divided into 
three sub-domains and a pulse was excited at the left boundary of the domain. The 
right boundary of the domain was terminated with PEC. A sample of the compu- 
tations is shown in Figure 7.6 where different frames of the propagating pulse, as 
well as the reflected one by the PEC, are shown. It seems that the coupling of 
the sub-domains was successfully implemented and there are no spurious reflections 
occurring at the interfaces. Here, a 17-node Chebychev grid was used in each sub- 
domain (including the fictitious points). The matching of the sub-domains (related 
to the determination of the field values at the fictitious points) at each interface was 
implemented such that it was dependent on the matching at other interfaces. The 
same experiment was also performed with the matching between two neighbor sub- 
domains (related to the determination of the field values at the fictitious points) not 
being dependent on the matching at other interfaces. The computations are shown 
in Figure 7.7 where different frames of the propagating pulse, as well as the reflected 
one by the PEC, are showm It is observed that both approaches of matching the 
interfaces give very similar results. Therefore, for the following experiments it is cho- 
sen to do the matching at each interface independently of the others. This makes the 
coding of a general sub-domain program much more simple and practical. 

The next simulation examines a 1-D domain divided into two sub-domains where 
one is free-space and the other one is PML. This is done to investigate how coupling 
between two sub-domains with different material properties is done. Following again 
the method proposed in [21], continuity of the fields as well as of their first time 
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Figure 7.6: Traveling pulse in 1-D domain terminated with PEC (matching at one 
interface depends on the others). 
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Figure 7.7: Traveling pulse in 1-D domain terminated with PEC (matching at 
interface is independent of the others). 
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derivative is enforced because the PML has the same permittivity e and permeabil- 
ity // with the free-space. However, in this case the continuity conditions are more 
complicated than the ones in the case of materials with no magnetic or electric con- 
ductivity. The reason is that in this case these conditions now involve both E and H 
fields coupled in the same equations. A sample of the simulation results are shown in 
Figure 7.8 where different frames of the propagating pulse are shown. As expected, 
when the pulse enters the PML it gets absorbed without creating any reflections. 



Figure 7.8: Traveling pulse in 1-D domain terminated with PML. 

Finally a general multi-domain code was written that can divide a domain into a 
number of sub-domains which have the same number of points. This code works only 
for homogeneous regions, i.e., all sub-domains must have the same material properties 
in addition to zero electric and magnetic conductivities. However, the code can be 
generalized for different types of materials. To illustrate the high accuracy of the 
Chebychev collocation method, a free-space domain was analyzed. The domain was 
two meters long and was excited at the left boundary with a Gaussian pulse that 
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had significant frequency components up to 1.5 GHz (10% bandwidth). On the right 
boundary the domain was terminated on PEC, but the simulation time was chosen 
so that the pulse does not reach this end creating reflections. This problem was 
analyzed first using the Chebychev collocation method for a different number of sub- 
domains and grid points, and second, using a second-order accurate FDTD scheme. 
The results are shown in Figure 7.9 and we observe that additional subdomains do not 
necessarily provide better accuracy. Also, the high accuracy of the spectral methods 
compared to the second-order FDTD scheme is again observed. 



Number of grid points 


Figure 7.9: Maximum error in the solution of the 1-D TM mode. 


V. Conclusions 

In this report, spectral methods have been applied for the first time in the context 
of electromagnetics. The 1-D Maxwell’s equations were successfully solved by these 
methods, and the PML was used to simulate open-space problems. The numerical 
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simulations illustrated that spectral methods outperform the standard FDTD algo- 
rithm in terms of accuracy, as they exhibit exponential convergence, i.e., the error 
drops exponentially as the number of discretization points increases. Furthermore, 
these methods exhibit negligible dispersion and therefore can simulate accurately 
electrically large domains. The main drawback of them is their restrictive stability 
criteria which require a very small time-step in order for them to be stable. However, 
there have been a few proposed methods that deal successfully with this problem and 
accomplish to relax the time-step restrictions. 

Moreover, the concepts of interior domain grids and fictitious points were intro- 
duced. These principles were applied in the multi-domain formulation that was pre- 
sented. A multi-domain approach has the advantage of using smaller differentiation 
matrices and it can also use larger time-step due to the fact that each sub-domain is 
formed by a number of points considerably smaller than the one of the entire domain 
method. This can lead to computational savings as a significantly larger time-step 
can be used without substantial loss of accuracy. 

The spectral methods seem to be very promising for the analysis of electromagnetic 
problems. They can handle efficiently large domains because of their high accuracy, 
exponential convergence, and negligible dispersion and dissipation. Our ultimate goal 
is to apply and extend these methods initially to generic 2-D and 3-D problems and 
then to radiation problems. 
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